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ABSTRACT 


The  restoration  of  images  and  the  enhancement  and 
detection  of  targets  in  cluttered  background  are  the  subjects 
of  this  research.  The  statistical  approach  is  used  in  order 
to  exploit  temporal  as  well  as  spatial  image  redundancies. 

The  images  are  modeled  as  a homogeneous  random  field. 

An  autocorrelation  function  and  a method  of  parameter 
identification  are  proposed.  Experiments  with  several 
pictures  are  presented  to  validate  the  model. 

An  analysis  of  two-dimensional  recursive  filters  is 
presented.  A three-dimensional  recursive  filter  is  developed 
which  exploits  the  spatial  as  well  as  the  temporal  image 
redundancies. 

A class  of  hybrid  filters  is  proposed  which  improves 
the  performance  of  the  recursive  filters.  Several  experi- 
ments with  pictures  are  presented  to  show  the  ability  of 
the  hybrid  filters  in  picture  restoration. 

A detector  is  developed  for  purposes  of  target  extrac- 
tion from  cluttered  background  images.  The  detection  is 
independent  of  the  target  shape. 

A simulation  of  the  target  detection  and  tracking 
problem  is  presented.  The  target  is  tracked  from  frame  to 
frame  by  means  of  a conventional  Kalman  filter,  which  uses 
the  image  filter  as  the  measurement  device. 
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I.  INTRODUCTION 

A.  PROBLEM  DEFINITION 

This  research  involves  two  aspects  of  digital  image 
processing.  One  is  the  restoration  of  images  degraded  by 
white  Gaussian  observation  noise.  The  other  is  the  enhance- 
ment and  detection  of  targets  immersed  in  cluttered  back- 
ground images.  The  statistical  approach  is  used  in  order 
to  exploit  temporal  as  well  as  spatial  image  redundancies. 
The  emphasis  is  directed  towards  the  design  of  recursive 
and  hybrid  Bayesian  filters. 

In  the  recent  past  considerable  attention  has  been 


devoted  to  the  application  of  Kalman  filtering  to  smoothing 
of  observation  noise  in  image  data.  Habibi  [4]  first  sug- 
gested a two-dimensional  recursive  image  estimator  as  an 
extension  of  the  one-dimensional  Kalman  filter.  Two  other 
similar  extensions  were  proposed  [3],  [5].  It  has  been 
shown  [10]  that  these  filters  do  not  preserve  the  optimality 
of  the  one-dimensional  Kalman  filter.  The  most  complete 
study  of  optimal  two-dimensional  Kalman  filtering  has  been 
performed  by  Woods  and  Radewan  [11] . They  point  out  that  the 
generalization  of  the  Kalman  filter  to  two  dimensions  can 
be  done  optimally  with  an  extremely  high  dimensional  state 
vector,  which  has  dimensions  on  the  order  of  MN  (M  = order 
of  the  filter,  N = width  of  the  image) . Panda  and  Kak  [9] 
succeeded  in  deriving  a vector  dynamic  model  that  generates 


i 


the  same  random  field  of  Habibi 1 s model.  Since  this  model 
is  recursive  in  only  one  index,  the  one-dimensional  Kalman 
filter  is  applied  and,  of  course,  optimality  is  preserved. 
Again,  as  indicated  by  Woods  and  Radewan  [11] , an  extremely 
high  dimensional  state  vector  is  used,  which  has  the  same 
dimension  as  the  width  of  the  image. 

Several  techniques  have  been  investigated  for  background 
clutter  suppression  and  target  enhancement.  Statistical 
non-recursive  spatial  filters  [6]  are  suggested  for  background 
clutter  suppression  and  enhancement  of  targets  of  known  shapes 
Two-dimensional  Kalman  filtering  [5]  is  also  proposed  for 
the  case  of  targets  with  arbitrary  shapes.  These  filters 
assume  that  the  image  gray  level  can  be  decomposed  into  three 
additive  components:  target,  background  and  white  observa- 
tion noise.  The  statistical  differences  of  these  three 
components  are  used  in  order  to  extract  the  target,  which 
is  the  desired  information. 

B.  RESEARCH  OBJECTIVES 

In  the  context  of  picture  restoration  and  target  detection 
as  defined  above,  several  specific  research  objectives  have 
been  identified. 

a)  Since  the  model  of  images  plays  a fundamental  role 
in  the  statistical  approach  of  image  processing,  several 
experiments  with  real  life  pictures  will  be  accomplished 
in  order  to  validate  existing  models  [1-2]  and  to  suggest 


new  ones . 


b)  Generalization  of  the  one-dimensional  Kalman  filter 
to  two  dimensions  results  in  excessive  computation  loads. 

On  the  other  hand,  the  filters  [3-5]  are  very  simple  and, 
although  it  is  known  that  they  are  not  optimal,  it  was  not 
determined  yet  just  how  far  they  are  from  optimality.  There- 
fore, these  filters  will  be  compared  against  the  optimum 
non-recursive  filter. 

c)  The  existing  recursive  image  filters  exploit  only  the 
spatial  correlation,  therefore  a three-dimensional  recursive 
filter  will  be  developed  in  order  to  take  advantage  of  the 
correlation  in  time. 

d)  The  possibility  of  improving  the  performance  of  the 
sub-optimum  recursive  filters  by  using  the  non-causal  (with 

a 

respect  to  the  direction  of  recursion)  observations  closest 
to  the  estimated  pixel  (picture  element)  will  be  investi- 
gated. The  result  is  a hybrid  filter  in  the  sense  of  using 
some  observations  recursively  and  others  non-recursively . 

e)  An  optimum  decision  rule  will  be  derived  for  purposes 


of  background  suppression  and  target  enhancement  and  subse- 
quent threshold  detection. 

f)  A conventional  Kalman  filter  will  be  constructed  to 
track  the  target  centroid  from  f rame-to-f rame . 


Second,  the  results  of  several  experiments  with  real  life 
pictures  are  presented.  Third,  a model  for  a picture  and  for 
pictures  sequenced  in  time  is  introduced. 

In  Chapter  III  the  sub-optimum  two-dimensional  recursive 
filters  [3-5]  are  analyzed.  The  exact  computation  of  the 
error  variance  is  particularly  important  to  those  filters, 
because  they  either  do  not  compute  the  error  variance  in 
the  calculation  of  the  gains,  or  use  an  approximation.  An 
exact  method  is  developed  to  compute  the  error  variance  of 
those  filters.  Using  such  a method,  the  filters  are  compared 
among  themselves  and  against  the  optimum  non-recursive 
interpolator,  constrained  to  the  same  data  set.  A recurs- 
ive filter  is  also  introduced  that  is  essentially  the 
same  as  [5] , but  the  computation  of  gains  is  accomplished 
without  approximation. 

In  Chapter  IV  a three-dimensional  recursive  filter  is 
developed.  This  filter  estimates  the  pixel  gray  level  of 
pictures  sequenced  in  time  by  using  recursively  the  obser- 
vations of  the  estimated  frame,  as  well  as  those  in  the  pre- 
vious frames.  It  is  an  extension  of  the  two-dimensional 
recursive  filter  [3].  Numerical  results  are  presented  to 
evaluate  the  improvement  resulting  from  the  exploitation  of 
the  correlation  in  time. 

Chapter  V introduces  a new  class  of  image  filters, 
called  hybrid  filters.  These  filters  are  smoothers  that 
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combine  optimally  the  estimate  of  the  recursive  filters 
(two  or  three  dimensions)  with  an  arbitrary  set  of  "future" 
observations . Theoretical  comparison  between  the  hybrid 
filter  and  the  recursive  [3]  and  non-recursive  [6]  filters 
is  accomplished.  Experiments  with  real  life  pictures 
are  also  presented  for  the  purpose  of  comparison  of  these 
filters . 

In  Chapter  VI  an  optimum  decision  rule  is  developed  to 
detect  targets  immersed  in  cluttered  background  images. 

The  target  is  considered  as  another  texture  statistically 
distinct  from  the  background  texture.  The  decision  is 
made  pixel  by  pixel  and,  therefore,  it  is  independent  of 
the  target  shape,  but  it  can  also  be  applied  to  targets  with 
known  shapes.  The  decision  is  based  on  the  observation  of 
the  pixel  gray  level  and  the  background  prediction  for  the 
pixel.  This  prediction  is  given  by  the  recursive  or  hybrid 
filters  and  may  also  be  given  by  a non-recursive  filter. 

Some  special  cases  are  worked  out  and  its  performance 
evaluated.  The  image  is  modeled  as  a weighted  addition  of 
three  components:  target,  background  and  observation  noise. 

In  Chapter  VII  a conventional  Kalman  filter  is  constructed 
to  track  the  target  centroid  (or  other  points)  from  frame 
to  frame.  The  target  dynamics  in  the  picture  is  modeled.  The 
detector  developed  in  Chapter  VI  feeds  the  tracking  filter 
with  the  observations  of  the  centroid  spatial  coordinates. 

The  observation  error  is  analyzed. 
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Chapter  VIII  presents  several  results  of  the  simulation 
of  a complete  target  detection  and  tracking  problem,  using 
computer  generated  images.  Both  recursive  and  hybrid 
filters  are  used  and  compared. 

The  final  chapter  summarizes  the  results  of  this 
investigation  and  presents  the  conclusions  and  suggestions 
for  further  study. 
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II.  IMAGE  MODELING 


A.  INTRODUCTION 

Our  objective  in  this  chapter  is  to  obtain  a statisti- 
cal model  for  pictures.  Knowing  part  of  a picture,  one 
can  generally  draw  certain  inferences  about  the  remainder; 
or,  knowing  a sequence  of  frames,  one  can,  on  the  average, 
make  a good  guess  or  prediction  about  the  next  frame.  From 
a statistical  viewpoint,  similarity  between  adjacent  pixels 
(picture  elements) , or  , frame- to- frame  similarity  represent 
a high  level  of  intraframe  or  interframe  correlation. 

Experimental  evicence  [ 1 ] — [ 2 ] indicates  that  a mono- 
chromatic image  can  be  modeled  by  specifying  its  value  (gray 
level  x(m,n)  at  each  spatial  coordinate  (m,n) . An  ensemble 
of  such  images  can  be  modeled  by  interpreting  x(m,n)  as  a 
random  field. 

In  this  chapter,  first,  we  will  address  the  mathematical 
problem  of  finding  the  dynamic  model  of  random  fields,  given 
the  autocorrelation  function.  Second,  we  will  introduce  a 
model  for  a picture  and  for  pictures  sequenced  in  time. 

Third,  we  will  present  experimental  results  in  order  to 
validate  the  proposed  model. 

B . DYNAMIC  MODEL 

Assume  the  mean  and  autocorrelation  function  of  a 
homogeneous  (wide  sense  stationary)  random  field  are  given. 
Since  we  are  assuming  knowledge  of  the  mean,  for  convenience, 


we  assume  the  random  field  has  zero  mean,  therefore  auto- 
correlation and  autocovariance  are  the  same. 

A highly  desirable  characteristic  of  any  dynamic  model 
is  to  have  it  excited  by  uncorrelated  noise.  The  impulse 
response  method  always  has  such  characteristics.  The  diffi- 
culty with  this  method  is  that  the  power  spectrum  density 
has  to  be  factored,  and  there  is  no  mathematical  theorem 
for  factorization  in  the  multi-dimensional  case.  Since 


there  exists  considerable  evidence,  in  refs.  [ 1 ] — [ 2 ] and 
in  this  research,  that  images  are  well  modeled  by  auto- 
correlation functions  with  separable  kernels,  the  difficulty 
of  factorization  will  not  arise  and,  therefore,  the  impulse 
response  method  will  be  used. 

The  impulse  response  method  is  shown  in  figure  2.1. 
Assume  there  exists  a system  transfer  function  H(z),  such 
that,  when  driven  by  white  noise  W(n),  the  resultant  output 
X(n)  is  a random  process  with  the  desired  autocorrelation 
function  R(k) . Using  input-output  relationships  between 
power  spectrum  densities  (PSD) : 

S ( z)  = S.  (z)H(z)H(z'1) 
o 1 

Since  the  input  is  white  noise,  its  PSD  is  constant, 


Figure  2 . 1 

Filter  Response  Method 


Therefore,  to  compute  H(z)  we  have  to  factor  the 
right  hand  side  of  the  previous  equation.  This  is  always 
possible,  in  the  one  dimensional  case,  provided  SQ(z)  is 
a ratio  of  polynomials. 

In  the  remainder  of  this  section  we  present  some  exam- 
ples to  illustrate  the  method,  as  well  as  for  later  use. 
Example  1 

Consider  the  autocorrelation  function: 

R(k)  = o2  e-a  ' k ' k = 0,±1,±2,  ... 

let 


R(k)  = a2  p|kl  (2.1) 

I 


Take  the  two-sided  Z-transform  of  R(k) , 


00 


z (R (JO  ) = l R(k)  z" 


k=-°° 


? 2 2 -12-2 
a (...+  p z + pz  + 1 + pz  + pz  + 


Z (R(k) ) = 


2 2 
o (1-P  ) 

( 1— pz  1) (1-pz) 


(2.2) 


Si(z)H(z)H(z“1) 


2 2 
a (1-P  ) 

(1-pz-1) (1-pz) 


To  have  a stable  model  the  poles  of  H(z)  must  be  inside 
the  unit  circle,  therefore  we  choose: 


1-pz 


Si(z)  = a2(l-p2) 


X(z) 

W(z) 


Thus,  the  dynamic  model  is: 


X(n)  = pX(n-l)  + V7(n) 


(2.3) 


Initializing  equation  (2.3)  with  X(0),  and  exciting 

2 2 

with  white  noise,  having  variance  a ( 1— p ) , a random  process 
X(n)  is  generated  having  the  desired  autocorrelation 
function. 
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Example  2 

Consider  the  autocorrelation  function: 


R(k)  = 7- — ytt—  ■T[(l-p.2)pJlc|+1 

pl~p2  (1+PiP2  2 1 


-(1-p,  ) P- 


1*1+1, 


(2.4) 


where : 


-a. 


-a. 


= e 


= e 


k = 0 , ±1 , ±2  , 


Using  results  of  example  1,  the  Z-transform  of  R(k)  is: 


Z(R(k))  = 


(1-p^z  1) (l-p2z  1) (l-pxz) (l-p2z) 


(2.5) 


where : 


a (l-p.p,)  n o 

q - -i+pJr-^^i  )a"p2 } 


Let's  choose: 


H(z)  = 


..  X(z) 


W(z) 


(1-P1z_1)  (1-P2z  1) 


Si(z)  = Q 


Thus,  the  dynamic  model  is: 
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X (n)  = (P1+P2)x(n-1)  - P1P2X(n-2)  +W(n)  (2.6) 

It  can  be  easily  shown  that  a particular  case  of  this 
model,  when  p^  approaches  p2,  is: 

R(k)  = a2  p (1 + 6 |k| ) (2.7) 


where 


6 = 


1 - p; 
1 + p' 


p = e 


-a 


and  the  dynamic  model  is: 


X (n)  = 2pX (n-1)  - p X(n-2)  + W(N) 


(2.8) 


Si(z) 


J2(l  - 


1 + P 


2 . 3 

Q^_) 


From  equation  (2.7)  we  can  see  that  this  autocorrela- 
tion function  has  zero  derivative  at  the  origin  (k  = 0) 
and  an  inflection  point  at  k * 1/a,  provided  that  p > 0.7 
or  6 z a. 

The  model  given  by  equation  (2.4)  is,  therefore,  a more 
general  case  and  includes,  as  particular  cases,  the  models 
given  by  equations  (2.1)  and  (2.7).  In  figure  2.2  some 
curves  of  equation  (2.4)  are  shown  for  the  same  "correlation 
time"  ( the  point  where  the  correlation  is  37%  of  the 


maximum) . 


AUTOCORRELATION  R(k) 


correlation  time 

(1)  o2  = ° 

(2)  px  = p2 

(3)  P]_  = 2p2 


SHIFT 


Figure  2.2 


Autocorrelation  Functions 


Example  3 

Consider  a two-dimensional  autocorrelation  function, 
separable  in  the  two  dimensions: 


R ( i ..  j ) 


(2.9) 


where 


av  -ah 

Py  — ® © , i , J - Or  —1  , • 


To  find  the  Z-transform  we  take  advantage  of  separa- 
bility and  use  equation  (2.2)  from  example  1. 


Z(R(i,j))  = 


a2(l-pv2)(l-ph2) 


-1 


-1, 


(l“Pvzi  ) (l-P„Zi) (1-Pk*,) 


h 2 


v‘1'  vx 


(2.10) 


Let's  choose; 


H(zlfz2) 


X(z1,z2) 


w(zx,z2) 


Si (zi'z2}  = a2(l-pv2) (1-Ph2) 


Thus,  the  dynamic  model  is: 


X (m,  n) 


P X (m-1  ,n)  +p,  X (m,n-l)  - p p.X  (m-1  ,n-l)  +W  (m,n) 


«•  -* 


Example  4 

Similar  to  the  previous  example,  consider  a separable 
autocorrelation  with  the  kernels  of  example  2: 


R(irj)  = a Kv(i)  Kh(j; 


(2.12) 


where : 


Kv(i) 


1 Tit  2 . I i i +1  2,  I i ! +1 , 

(p.  -p.  ) (1+p.  p„„) 1 (1-p2v  ) plv  ■{1"piv  )p2v  ] 


lv  2v  ' lv  2v 


Kh  (j ) = same  equation  with  i -*■  j 

v h 


The  Z-transform  is  given  by: 


00  00 

Z (R(i, j ) ) = 11  R(i»j)  z^1  z2_j 

— OO  j =— CO 


00  00 


a l [ l Kh( j) z2  J]Kv(k)z1 
i = -00  — 00 


a zn  (Kv(i) ) z2 (Kh( j) ) 


-i 


Using  equation  (2.5)  it  can  be  easily  verified  that  the 
transfer  function  is: 
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(l“plv2l  ) (l-p2vZl  )(1-plhz2  )(1-p2hZ2  5 


(z^ , z2)  - o Qv 


where : 


v 


(1_Plvp2v)  „ _ 2W1  . 2 , 

T-— — -(1-Piv  ) d-P2v  ) 


lvy2v 


= same  equation  with  v -►  h 


Thus,  the  dynamic  model  is: 


X(m,n)  = (plv+p2v)X(m-l,n)  - p1vp2vX (m-2 ,n) 


+ (Plh+P2h) X (m,n-l)  - plhp2hX(m,n-2) 


(piv+p2v)(pih+p2h^  X ,n-1^ 


+ plvp2v(olh+p2h)X(m"2'n"1) 


+ plhp2h(plh+p2v)X(m“1'n~2) 


- P1vP2vPlhP2hX(m-2,n-2)  + W(m,n) 


(2.13) 


(2.14) 


Observe  that  example  3 is  a particular  case  of  this, 


where  p2v  = p2h  = 0, 


2< 


■HI 


sc** 


To  simplify  notation,  from  now  on,  we  will  adopt  the 
following  convention: 


X(m,n)  = A X + W(m,n) 


where: 


A = row  vector  of  coefficients 


X = column  vector  of  adjacents 


In  the  case  of  example  3: 


A = 


[p  p.  -p  p,  ] 

v h v h 


X = 


X (m-1  ,n) 

X (m,n-l) 

X (m-1 , n-1) 


(2.15) 


Example  5 

Consider  a three-dimensional  separable  autocorrelation 
function  having  kernels  like  examples  1 and  2: 


X 


(l-plvzi  ) (1-p2vZl  * (1-pIhZ2  * (1-p2hZ2  * {1-ptZ3 


(2.17) 


Silzl'z2'z2)  = °2  QvQh  (1~pt2) 

where  Qv  and  are  the  same  as  in  example  4. 

Thus,  the  dynamic  model  is: 

X (m, n, t)  = A X + W(m,n,t)  (2.18) 


where : 

A is  a row  vector  whose  elements  are  the  coefficients 
of  eq.  2.14,  followed  by  these  same  coefficients 
multiplied  by  pfc,  and  the  final  element  is  pfc. 

X is  a column  vector  whose  elements  are  the  adjacent 
pixels  of  eq.  2.14,  followed  by  these  same  pixels 
located  on  the  previous  frame,  and  the  final  element 
is  X (m,n, t-1) . 


In  figure  2.3  the  adjacent  pixels  are  shown. 
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Figure  2.3 

Adjacent  pixels  of  example  5 

C.  AN  AUTO- REGRESSIVE  MODEL 

In  this  section  we  present  a stochastic  model  for 
images.  The  basic  assumption  is  that  images  constitute  a 
homogeneous  random  field  with  zero-mean  (or  known  mean) 
and  known  autocorrelation  function,  separable  in  the  inde- 
pendent dimensions.  These  assumptions  will  be  validated 
by  experiments  in  the  next  section. 

1 . Autocorrelation  Function 

The  autocorrelation  function  chosen  is  intended  to  be 
general  enough  to  include  most  of  the  models  used  in  other 
research  (2]  - [9]  , as  well  as  to  best  fit  the  experimental 
functions  measured  in  this  research  and  [1]— [2] . The 
hypothesis  of  separability  allows  us  to  examine  the  auto- 
correlation kernel  by  kernel.  Assume  that  the  process  is 


modeled  by  a kernel  of  first  order,  as  in  Example  1. 
The  dynamic  model  is 


X(n)  = p^X(n-l)  + W^(n) 


If  this  model  is  valid,  the  modeling  error  W-^(n)  must  be 
uncorrelated  noise.  This  can  be  easily  verified  by  computing 
the  autocorrelation  of  W^Cn)  as  follows.  We  compute  W^(n) 
at  each  point  by 


Wx(n)  = X (n)  - p^X (n-1) 


and  then  compute  the  autocorrelation  of  this  sequence. 
Computing  the  autocorrelation  of  W^(n)  is  a very  good 
test.  Assume  that  W^(n)  turns  out  to  be  correlated. 

In  this  case  we  can  model  W^(n)  again  by  a first  order 
kernel : 


W^(n)  = p2W1(n-l)  + W2 (n) 


Assume  that,  after  repeating  the  same  measurements  for 
W2 (n) , we  conclude  that  W2(n)  is  uncorrelated  noise. 

In  figure  2.4  it  is  shown,  in  the  Z-domain,  the  operations 
that  we  have  performed. 

Thus,  the  overall  transfer  function  is : 


» 
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w1(z) 

H 2(z) 

■■ 

Figure  2.4 
Second-order  model 


(1-P1z  ) (1-P2Z  ) 


The  second-order  dynamic  model  is: 


! (n)  = (P1+P2) X(n-l)  - p1P2X(n-2)  + W2(n) 


(2.19) 


We  have  now  a model  driven  by  white  noise,  which  is  a 
second-order  difference  equation. 

The  process  generated  by  this  model  has  an  autocorrela- 
tion function  given  by  equation  (2.4),  as  we  will  demonstrate 

E(W  2(n))  = Q- 


Z(R(k) ) 


Q2  H ( z)  H (z  ■l) 


where : 


since 


(1-P1z  1) (l-p1z) (1-P2z  l)  (1-P2z) 

AQ2 + BQ2 

(1-p^-1)  (l-p1z)  (l-p^'1)  (l-p2z) 

P1 

(P1“P2)  (l-P-^) 

~p2 

(p1-p2)  ( 1— P ^P  2) 


z'H 


(1-pz-1) (1-pz) 


v _ 1 Jkj 


1-P‘ 


R(k)  = Q0[ 


2 1 , 2 P1 
l“Pi 


1-P 


B ,1*1] 


2 2 


Using  R(0)  = a (variance  of  X(n)) 


2 2 2 ^ P l P 9 ) 

o2  = «2(1-Pl2)(i-P22)  (1+p^) 


(2.20) 


R (k)  = 


2 . _ k 


(P1-P2) (1+p1p2) [pl(1"p2  )pl'  “p2(1“pl  )p2 


(2.21) 


where : 


A1  = p^2  (p2*p3)  (-p2p3)  (l-p22)  (l-p32) 

4 

A2  = -P22(P1“P3)  (1-P103)  (1-D12)  (1-P32) 

A3  = P32  (prp2}  (1"plp2}  (1_0l2)  {1_P22) 


X(n) 


The  general  expressions  for  a k-order  model  are: 
Dynamic  model: 


~ ^pl+p2+*  ’ * pk^  X ~ ( plp2+plp3+*  • * pk-lpk)  X (n-2 ) 

+ (P]_P2P3+- * * Pk-2Pk-lPk^  X^n-3) 
k+1 

+ (-1)*  ■lp1p2.  . .PjcX(n-k)  + Wk(n)  (2.25) 

Autocorrelation  function: 


R(i)  = 


2 Alpl 


i| 


+ A2p2 


Vk1 


A]_  + A2  + 


+ A, 


(2.26) 


Variance  of  Wk(n): 


(2.27) 


The  expressions  for  the  A's  are  easily  inferred 
from  equation  (2.23). 

In  example  3 an  auto-regressive  model  is  presented 
for  a picture  with  an  autocorrelation  whose  kernels 
in  the  two  spatial  dimensions  are  the  first  order  versions 
(k=l)  of  equation  (2.26).  In  fact,  that  model  is  used  in 
most  research  [21  - [9]  with  images.  In  example  4 
an  image  model  is  presented  whose  kernels  are  the  second- 
order  versions  (k=2)  of  equation  (2.26).  The  order  of  the 


kernels  doesn't  have  to  be  the  same.  Depending  on  the 
particular  picture  to  be  modeled,  the  order  is  chosen  for 
each  direction. 


In  example  5 an  auto-regressive  model  is  presented 
for  pictures  sequenced  in  time.  In  that  case  we  chose 
second-order  kernels  in  the  spatial  directions  and  first- 
order  in  the  time  direction. 

At  this  point,  some  special  features  of  the  proposed 
autocorrelation  functions  given  by  equation  (2.26)  should 
be  emphasized: 

a)  It  includes  as  a particular  case  the  first 
order  model  which  is  used  in  most  research  [2]  - [9]. 

b)  It  fits  much  better  many  situations  where  the 
first  order  model  is  a poor  approximation. 

c)  The  resultant  dynamic  model  is  auto-regressive 
and  driven  by  white  noise. 

d)  The  cascade  feature  of  whitening  the  modeling 
error  is  quite  simple  and  adequate  for  on  line  parameter 
identification.  This  feature  will  be  exploited  in  the  next 
item. 

2 . Parameter  Identification 

Given  that  the  autocorrelation  function  of  an 
image,  or  a time-frame  of  images,  is  assumed  to  have 
I separable  kernels  as  in  equation  (2.26),  the  next  problem 

| is  to  find  its  parameters.  Let's  discuss  the  modeling  of 

just  a picture,  since  the  extension  for  the  time-frame  is 
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straightforward.  Having  a specific  picture,  we  want  to 
identify  the  parameters  of  the  model.  To  compute  statistics 
from  just  one  realization  of  the  process  we  must  assume 
the  process  is  ergodic,  therefore  we  consider  the  picture 
as  a typical  realization  and  the  averaging  in  space  is 
equivalent  to  statistical  averaging.  Also  the  size  of  the 
picture  has  to  be  sufficiently  large. 

Those  conditions  are  quite  severe,  and  we  realize 
that  the  process  has  to  be  restricted  to  a particular 
class  of  pictures  with  common  properties.  Another  way 
of  looking  at  this  is  to  visualize  the  picture  as  an  array, 
for  example,  of  100  x 100,  therefore  with  10,000  pixels. 

Each  pixel  can  have,  say,  100  distinguishable  brightness 
values.  The  number  of  different  pictures  is  10010'000. 

With  such  enormous  numbers  of  pictures,  even  if  we  were 
able  to  compute  the  mean  and  autocorrelation,  such  moments 
would  not  be  enough  information.  It  would  be  necessary 
to  know  statistics  of  enormously  higher  orders.  The  atten- 
tion, therefore,  has  to  be  focused  on  local  redundancy, 
encompassing  only  a few  pictures  closely  related.  To  have 
a model  as  accurate  as  possible,  based  only  in  the  first 
and  second  moments,  we  will  need  some  degree  of  adaptation. 

Our  method  is  based  on  the  cascade  feature  of 
equation  (2.26).  Let's  explain  the  method  in  detail. 

Assume  a picture  is  given  with  negligible  measurement 
noise.  We  proceed  in  the  following  steps. 


Step  1 - Compute  the  mean  and  remove  it  from  the  picture. 
Step  2 - Compute  the  autocorrelation  in  the  horizontal 
direction,  with  a number  of  displacements  equal  to  20%  of 
the  number  of  columns.  There  is  no  use  to  go  much  further 
than  this,  because  the  accuracy  in  the  averaging  is  reduced. 
Repeat  the  same  procedure  in  the  vertical  direction.  This 
step  can  be  accomplished  using  the  Fast  Fourier  Transform 
algorithm  (FFT) . 

Step  3 - Having  the  measured  autocorrelation  in  the  hori- 
zontal and  vertical  directions,  assume  that  equation  (2.26) 
is  adequate  and  gives  a nice  fit.  Start  with  a first-order 
kernel  and  compute  and  such  that  the  mean-squared 

error  in  the  fitting  is  minimized.  The  result  is  a model 
like  example  3: 

X(m,n)  = X + W^(m,n) 

where : 


-1  [plv  plh  "plvplh] 

XT  = [X (m-1 ,n)  X (m,n-l)  X(m-l,n-l)] 


Step  4 - At  this  point  we  have  a model  for  the  picture  that 
is  first-order  in  both  spatial  directions.  To  test  the 
goodness  of  this  model  the  obvious  figure  of  merit  is  the 
modeling  error  W^(m,n).  If  the  model  is  adequate,  W^(m,n) 


■ 
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should  be  uncorrelated  noise.  To  test  this,  compute  the 
autocorrelation  of  W^(m,n)  as  in  step  2.  If  it  is  corre- 
lated in  both  directions,  it  means  the  first-order  model  is 
not  good  enough  and  we  have  to  proceed  to  a second-order  model . 
It  could  happen  that  we  are  satisfied  with  one  direction, 
but  not  with  the  other,  in  this  case  we  have  to  proceed 
to  a second-order  only  in  such  directions.  Assume  here 
that  W^trrwn)  is  correlated  in  both  directions. 

Step  5 - Repeat  step  3 for  W^(m,n)  in  order  to  find  p 2v 
and  P2h*  The  model  for  W^(m,n)  is: 

W1(m,n)  = A2  W1  + W2(m,n) 


where : 


—2  * [p2v  p2h  “p2vp2h] 

= [W^(m-l,n)  W^(m,n-1)  W^(m-l,n-l)] 

The  second-order  model  for  the  picture  is  given  by 
equation  (2.14)  in  example  4. 

Step  6 - Repeat  step  4 for  W2(m,n).  If  it  is  sufficiently 
uncorrelated,  the  procedure  stops  and  we  have  identified 
the  order  of  the  model  as  well  as  the  parameters.  Other- 
wise we  proceed  seeking  a higher  order  model. 
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NOTE:  Step  4 is  responsible  for  the  decision  to  proceed 

or  not  and  also  in  which  directions,  whether  horizontal, 
vertical  or  both.  The  decision  to  proceed  or  not  can  be 
implemented  only  based  on  the  variance  of  W(m,n),  by  observing 
its  reduction  after  the  second-order  modeling.  If  such  reduc- 
tion is  null  or  negligible,  it  means  the  original  model 
is  good  and  it  is  not  necessary  to  increase  the  order.  The 
decision  about  the  direction  to  be  improved  can  be  imple- 
mented by  observing  the  autocorrelation  of  W(m,n)  at  a few 
points,  say,  the  first  5 displacements. 

The  computation  of  the  coefficients  is  quite  simple 
because  we  have  reduced  the  problem  to  fitting  exponentials 
at  each  step.  Such  computation  is  accomplished  as  follows. 
Assume  the  measured  autocorrelation  in  one  direction  is 
given  by  Y (n) , n = 0,1,... K and  it  is  normalized  such  that 
Y(0)  - 1.  We  want  to  fit  an  exponential  minimizing  the 
mean-squared  error,  R(n)  = e an  = pn.  Taking  the  natural 
logarithm  of  R(n)  and  Y(n),  we  want  to  minimize: 

K 

l [In  R (n)  - In  Y(n) ] 2 
n=l 

K 

l [-an  - In  Y(n) ] 2 
n=l 
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Take  the  derivative  with  respect  to  a and  equate 


1 

-• 

( 


to  zero,  the  result  being: 

K 

[ n In  Y(n) 

n=l 

K 

r 2 

l n 

n=l 


K 

-6  l n In  Y (n) 

n=l 

K (K+l) (2K+1) 


and 


In  figure  2.5  a flowchart  to  implement  the  proposed 
method  of  modeling  is  shown. 

To  model  pictures  sequenced  in  time  the  procedure 
is  quite  similar.  In  this  case  we  have  a sequence  of  pic- 
tures, sequenced  in  time,  and  apply  the  method  adding  the 
extra  dimension  (time) . Example  5 presents  a model  that  is 
first-order  in  time  and  second-order  in  both  spatial 
dimensions . 

3 . Modeling  Of  Noisy  Images 

In  the  last  item,  a method  for  picture  modeling  was 
presented  in  which  the  picture  was  assumed  to  be  noise  free. 
Such  an  assumption  restricts  considerably  real  life  appli- 
cations where  there  is  no  knowledge  of  the  original  picture . 
This  restriction  can  be  removed  if  the  autocorrelation 
function  of  the  measurement  noise  is  known. 


1 


i 

* 


Li 


38 


i 


I 

• ~m 

In  what  follows,  we  present  the  steps  to  be  modi- 
fied in  the  previous  modeling  method,  in  order  to  take  into 
account  additive  zero-mean  white  measurement  noise,  with 
known  variance. 


The  measurement  equation  is : 
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Figure  2.6 

Autocorrelation  of  y(m,n) 


j 


X(m,n)  = X + W^(m,n) 

7(111,0)  = X (m,n)  + v(m,n) 

Substituting  (2.28)  into  (2.29): 

y(m,n)  = Ax  Y + v(m,n)  + W^nwn)  - Aj.  V 

where : 

yT  = [y(m-l,n)  y(m,n-l)  y(m-l,n-l)] 

VT  = ( v (m-1 , n)  v (m, n-1)  v(m-l,n-l)l 
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(2.28) 


(2.29) 


y(m,n)  = Y + u^nwn)  (2.31) 


Since  y(m,n)  and  are  known,  u^(m,n)  can  be 
computed,  and,  therefore,  its  autocorrelation  in  both 
directions. 

Let's  examine  the  autocorrelation  of  u. (m,n)  in 
the  horizontal  direction. 

! 

I 

E [u^ (m,n) u (m,n+j ) ’ = E[v^(m,n)v^(m,n+j) ] + E [W^ (m,n) (m,n+j ) ] 

E [ (m, n) (m, n+ j ) ] = E [v (m,n) v (m, n+ j ) ] - A^E [ v (m,n) v (m,n+ j ) ] 

- A^Etvfn^n+j)  v(m,n)  ] 

+ A^E [v (m, n) vT (m, n+ j ) ] A^T 

where: 


E [v (m,n) v (m,n+j ) ] = r 5 (0  , j ) 


Figure  2 . 7 

Autocorrelation  of  u^(m,n) 

The  procedure  above  is  easily  extended  to  calculate 
the  autocorrelation  of  W2(m,n),  as  well  as  higher  order 
modeling  errors. 
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u^nwn)  = W^nwn)  + v1(m,n)  (2.33) 

* 

Substitute  (2.32)  into  (2.33): 

ux (m,n)  = A2  + u2(m,n)  (2.34) 

where : 

u2(m,n)  = W2(m,n)  + v2(m,n)  (2.35) 

v2(m,n)  = v1(m,n)  - A2  ^ 

From  equation  (2.34)  we  can  calculate  u2(m,n),  " 

since  u1(m,n)  was  computed  by  equation  (2.31),  and  the 

coefficients  A 2 are  given  by  the  fitting  of  (i,j). 

To  continue  the  modeling,  the  autocorrelation 

of  W_(m,n)  is  derived  from  R (i,j)  using  equation  (2.35). 

U2U2 

D.  EXPERIMENTAL  RESULTS 

In  this  section  we  present  some  relevant  results  obtained 
with  real  life  pictures.  Several  pictures  were  analyzed 
and  divided  in  two  broad  classes: 

Picture  A - represents  a class  of  pictures  with  few 

I 

I details  (i.e.  a low  content  of  high  spatial  frequency 

| structure).  In  figure  2.3  a sample  of  such  pictures  in 

a three-dimensional  plot  is  shown,  where  the  height  (z-axis) 

is  the  gray  level.  This  kind  of  plot  is  very  useful  to 

i 

emphasize  small  details. 

< II 
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Picture  B - represents  a class  of  pictures  with  many 


i 


details.  Figure  2.14  shows  a sample  of  such  pictures. 

1.  Autocorrelation  Function 

The  proposed  model  assumes  that  pictures  have 
separable  autocorrelation  functions  whose  kernels  are  given 
by  equation  (2.26).  To  validate  this  assumption,  the 
autocorrelation  of  pictures  A and  B were  measured  and 
compared  with  that  of  the  proposed  model.  Figure  2.9  is 
the  measured  autocorrelation  of  picture  A,  and  figure  2.15 
is  the  same  for  picture  B.  Using  minimum  mean-squared 
error  criteria,  a first-order  model  was  fitted  for  both 
pictures,  as  shown  in  figures  2.10  and  2.16.  Such  a model 
didn't  work  well  for  picture  A,  as  can  be  seen,  comparing 
figures  2.9  and  2.10.  For  picture  B it  worked  quite  well, 
as  can  be  seen  comparing  figures  2.15  and  2.16.  Figure  2.11 
is  the  second-order  model  for  picture  A;  observe  that  it  is 
> a very  good  fit.  Based  on  these  results,  we  conclude: 

a)  The  proposed  autocorrelation  is  adequate  for 
picture  modeling. 

b)  The  first-order  model  used  in  other  research 
does  not  fit  well  pictures  of  class  A. 

c)  Pictures  with  few  details  are  best  fitted  by 
a second-order  model. 

d)  Pictures  with  many  details  are  adequately 
fitted  by  first-order  model. 


I 

I 

I 
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2 . Parameter  Identification 

The  proposed  method  of  parameter  identification 
was  applied  to  both  pictures.  Figures  2.12  and  2 J. 7 are  the 
autocorrelation  of  W1 (m,n)  for  pictures  A and  B,  respectively. 
The  autocorrelation  of  the  modeling  error  W1(m,n)  is  the 
measurement  of  the  "goodness "-of  the  first-order  model.  As 
predicted  before,  picture  A has  a poor  first-order  model,  since 
W^(m,n)  is  quite  correlated,  instead  of  white  noise,  as 
can  be  seen  in  figure  2.12.  On  the  other  hand,  picture  B 
has  a nice  first-order  model,  since  W^(m,n)  is  almost  white 
noise,  as  can  be  seen  in  figure  2.17. 

Applying  the  cascade  method  of  modeling  for  picture 
A,  the  second-order  modeling  error  W2(m,n)  was  calculated 
and  its  autocorrelation  plotted  in  figure  2.13.  Now  the 
modeling  error  W2(m,n)  is  quite  close  to  white  noise,  and 
the  second-order  model  is  adequate  for  picture  A,  as 
predicted  before. 

3 . Bandwidth  Compression 

In  order  to  present  another  piece  of  evidence  to  validate 
the  model,  it  was  applied  to  bandwidth  compression.  Picture  B 
is  quantized  in  256  gray  levels,  thus,  8 bits  are  required 
to  transmit  the  gray  level  of  each  pixel.  Exploiting 
redundancies  in  the  picture  and/or  in  a time  frame  of 
pictures  (like  television) , it  seems  possible  to  reduce 
the  number  of  bits/pixels  to  be  transmitted,  therefore, 
reducing  the  required  bandwidth. 


I 
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The  model  for  picture  B is: 


X(m,n)  = A X + W^(m,n) 

The  measured  variances  of  X(m,n)  and  W1(m,n)  are: 

VAR  (X)  = 49 

VAR  (WL)  = 1 

Also,  it  was  seen  that  W^(m,n)  is  almost  uncorre- 
lated noise.  Therefore,  instead  of  transmitting  the  picture 
X,  we  can  transmit  : 

W^(m,n)  = X(m,n)  -AX 

Since  is  close  to  white  noise,  and  has  much 
smaller  variance  than  X,  the  quantization  levels  for  W may 
be  much  smaller  than  for  X.  This  method  of  bandwidth 
reduction  is  called  DPCM  (Differential  Pulse-code  Modulation) . 

Of  course,  the  effectiveness  of  this  method  is 
strongly  dependent  on  the  "goodness"  of  the  model. 

Figure  2.18  presents  the  histogram  of  the  modeling 
error  W^.  It  can  be  seen  that  most  "energy"  of  is 
concentrated  around  zero.  In  figure  2.19  the  reconstruction 
of  picture  B is  presented,  where  was  transmitted  with 
only  2 bits/pixel  (average) , instead  of  the  8 bits/pixel 
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required  for  X.  Therefore,  even  with  a bandwidth  reduction 
of  4 times,  the  reconstructed  picture  yet  carries  most  of 
the  information  of  the  original. 
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III.  TWO-DIMENSIONAL  RECURSIVE  FILTERS 


A.  INTRODUCTION 

A wide-sense  Markov  (WSM)  sequence,  corrupted  by  addi- 
tive white  noise,  can  be  optimally  (in  the  linear  least 
squares  sense)  estimated  by  a recursive  estimator.  The 
fact  that  image  intensity  can  be  modeled  by  a generaliza- 
tion of  a WSM  sequence  (Chapter  II)  motivates  one  to  search 
for  a two-dimensional  Kalman  estimator  to  smooth  additive 
white  observation  noise  in  images.  In  references  [3] , [4] 
and  [5]  three  filters  are  proposed.  These  filters  use  the 
dynamic  model  in  equation  (3.2),  that  is  a particular  case 
(first  order  WSM)  of  the  more  general  model  proposed  in 
Chapter  II.  The  structure  of  these  filters  is  also  the 
same  (see  equation  3.4),  the  difference  being  the  method  of 
gain  computation.  The  reason  of  the  existence  of  three 
Bayesian  filters  with  the  same  structure,  but  different 
gains,  is  that  they  are  all  sub-optimum.  It  has  been 
determined  [10]  that. the  optimality  of  the  Kalman  filter 
is  not  preserved  when  generalized  to  two  dimensions,  with 
such  dynamic  models.  In  reference  [9]  is  proposed 
a vector  model,  described  by  a first-order  linear  n-dimensional 
vector  difference  equation,  that  is  recursive  in  one  parameter 
and  generates  the  same  random  field  used  in  [3-5] . With 
such  a model  the  Kalman  filter  is  applied  preserving,  of 
course,  optimality.  Unfortunately,  the  complexity  of 
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implementation  of  such  a filter  is  enormous,  due  to  the 
size  of  the  state  vector,  that  is  the  number  of  pixels  in 
one  row,  say,  256  or  even  100.  Also,  its  complexity  increases 
drastically  if  we  need  to  extend  it  to  three-dimensional 
random  fields,  as  in  the  case  of  time-frame,  which  will  be 
addressed  in  chapter  IV.  Similarly,  the  extension  to  higher 
order  WSM  random  fields  increases  its  complexity. 

Because  of  the  above  considerations,  we  will  examine 
more  carefully  the  sub-optimum  filters  in  [3-5],  due  to  their 
simplicity  and  adequacy  to  accommodate  an  extra  dimen- 
sion (time)  or  a higher  order  WSM  random  fields. 

In  this  chapter  we  will  compare  the  three  filters 
[3-5]  among  themselves,  as  well  as  against  the  optimum 
non-recursive  interpolator,  constrained  to  the  same  data 
set.  We  will  also  introduce  a recursive  filter  that  is 
essentially  the  same  as  [5],  but  the  computation  of  gains 
is  accomplished  without  approximation.  The  figure  of 
merit  will  be  the  error  variance,  since  this  is  the  cost 
function  used  by  all  three  filters.  Observe  that  filters 
[4]  and  [5]  compute  error  variance  in  order  to  calculate 
the  gains,  but  there  is  an  approximation  in  the  recursive 
equation,  and  this  may  result  in  big  errors  in  the  gains, 
as  well  as  in  the  computed  error  variance. 

In  the  remainder  of  this  chapter,  we  will  develop  an 
algorithm  to  compute  the  variance  of  the  estimation  error, 
in  order  to  compare  the  filters  having  the  structure  of 
equation  (3.4)  and  dynamic  model  of  equation  (3.2).  Since 
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filters  [4]  and  [5]  have  similar  algorithms  for  computation 
of  gains,  and,  because  filter  [4]  is  clearly  inferior  to 
[5],  only  the  former  will  be  used  in  the  comparison. 

B.  ALGORITHM  FOR  COMPUTATION  OF  ERROR  VARIANCE 

Given  the  dynamic  model,  the  filter  structure  and  the 
gain,  we  are  going  to  develop  an  algorithm  to  compute  the 
variance  of  the  estimation  error. 

1.  Dynamic  Model 

The  picture  is  modeled  by  a zero-mean  random  field, 
homogeneous,  with  autocorrelation  function: 


R(i, j ) 


(3.1) 


The  dynamic  model  is: 


X(m,n)  = pvX (m-l,n) +phX (m,n-l) -pvphX (m-l,n-l) +W (m,n) 

(3.2) 

where  the  gray  level  X(m,n)  has  the  autocorrelation  of 
equation  (3.1)  and  W(m,n)  is  the  random  forcing  input  (or 
modeling  error).  W(m,n)  is  zero-mean  white  noise,  uncorre- 
lated with  X (p,q) , for  all  pixels  (p,q)  in  region  Xm  n 
defined  in  figure  3.1. 

E (W2 (m, n) ) = a2(l-pv2) (l-ph2)  = Q (3.3) 
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2.  Filter  Structure 


The  filters  [3,5]  use  equation: 

A A A 

X(m,n)  = [1-K (m,n) ] [pvX (m-l,n)  + phX(m,n-l) 

A 

- PvPhX(m-l,n-l) ] + K(m,n)y (m,n)  (3.4 

where  y(m,n)  is  the  measurement  of  X(m,n): 

y(m,n)  = X(m,n)  + v(m,n)  (3.5 

where  v(m,n)  is  zero-mean  additive  white  measurement  noise 
uncorrelated  with  X(m,n),  and  having  variance: 

E ( v2 (m,n) ) = r 

A 

The  estimate  X(m,n)  can  be  seen  as  the  optimal 

A 

linear  combination  of  the  prediction  X (m,n) , and  the 

hr 

measurement  y (m,n) : 

-A  A 

X(m,n)  = [1-K(m,n) ]X  (m,n)  + K(m,n)y(m,n)  (3.6 

r' 

where : 

A A A A 

X (m,n)  = p X(m-l,n)  + o,X(m,n-l)  - p p,X(m-l,n-l) 


The  prediction  X (m,n)  carries  all  the  information 
about  X(m,n),  contained  in  the  "past"  measurements.  The 


gain  K(m,n)  is  computed  such  that  all  the  information  about 

A 

X(m,n) , contained  in  Xp  (or  the  "past"  measurements)  and 
y(m,n),  is  effectively  used.  If  that  happens,  the  estima- 
tion error 

A 

e(m,n)  = X(m,n)  - X(m,n) 

is  uncorrelated  with  the  measurements  used  to  estimate 
X (m, n)  . 

It  is  shown  [10]  that  the  estimation  error  can  not 
be  orthogonal  to  all  the  "past"  measurements,  but  may  be 
only  to  the  region  defined  in  figure  3.2,  which 

excludes  row  m and  column  n. 

As  a matter  of  fact,  filters  [3-5]  do  not  have 
estimation  error  orthogonal  to  all  measurements  in  Y , . 

although  incorrectly  stated  in  [3]  and  [4].  A simple  way  to 
prove  this  above  is  using  a counterexample.  Let's  examine 
the  estimation  of  X(2,3).  According  to  [3]  and  [4],  the 
estimation  error  of  this  pixel  should  be  orthogonal  to 
y(l,2).  Also  the  error  in  the  estimation  of  X(2,2)  can't 
be  orthogonal  to  y ( 1 , 2 ) , according  to  [10].  We  are  going 
to  verify  that  these  statements  are  incompatible. 


The  error  in  the  estimation  of  x(2,3)  is 


e(2,3)  = X (2 , 3 ) - X (2 , 3) 

Using  equations  (3.2),  (3.4)  and  the  fact  that 

W(2,3)  and  v(2,3)  are  uncorrelated  with  y(l,2),  it  can  be 
verified  that 


e(2,3)y(l,2)  = [1-K (2 , 3) [pve (1 , 3 ) y (1 , 2 ) 


+ Phe(2,2")y  (1,2)  - pvphe (l, 2) y (1, 2)  1 


Since  the  filters  reduce  to  the  one-dimensional 
Kalman  filter  for  the  first  row: 


e(l,3)y(l,2)  = e(l,2)y(l,2)  = 0 


Thus 


e (2 , 3)  y (1, 2)  = [1-K  ( 2 , 3)  ] p.  e ( 2 , 2)  y (1 , 2) 


* 


Therefore,  since  in  general  K(2,3)  ^ 1 and  e(2,2)  can  not 
be  orthogonal  to  y(l,2),  we  conclude  that  e(2,3)  is  not 
orthogonal  to  y(l,2). 

3 . Variance  of  the  Estimation  Error 

In  the  following  derivation  we  will  simplify 
notation  by  dropping  the  argument  (m,n) , where  this  does 
not  cause  confusion. 

The  variance  of  the  estimation  error  is  given  by: 

P(m,n)  = E [ (X-X)  2 ] = E(X2)  + E(X2)  - 2E(XX) 
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Let's  call: 


2 

E (X  ; = 

/\ 

E (XX)  = b 

P(m,n)  = + a ^ - 2b  (3.8) 

e 

a.  Variance  of  the  Estimate 

2 

To  compute  ag  , observe  that  X can  be  written 

as  a linear  combination  of  all  measurements  in  region 

Y (see  figure  3.2). 
m,n 

X (m,n)  = bi;Ly(l,l)  + ...  blqy(l,q)  + ...  blnY(l,n)  + 
bplY(P,1)  + bpqY(p,q)  + bpnY(p,n)  + 


Substituting  in  (3.9): 


X(m,n) 


T 

= BaY 


The  estimate  variance  is: 

2 T T 

a0  = E(X  ) = B E(YY  )B 


(3.10) 


or 


a 

e 


2 


T 

B Ry  B 


(3.11) 


where : 


Ry  = E(YYT) 

Using  equations  (3.1)  and  (3.5),  it  can  be 
easily  seen  that  the  elements  of  the  autocorrelation  matrix 
Ry  are  given  by: 


Ry ( i » j ) 


R (k-p , Z-q)  (k,£)  ^ (p,g) 

< 

a2  + r (k,£)  = (p,q) 


(3.12) 


where  i is  the  sequential  number  for  the  pairs  (k,i)  and 
j for  (p , q)  as  follows 
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4 . Calculation  of  the  Coefficients 

To  find  the  vector  of  coefficients  B,  observe  in 
equation  (3.9)  that  when: 


y(ifj)  = { 


(i»j)  = (p,q) 


(i » j ) ? (p,q) 


(3.15) 


The  result  is: 


X(m,n)  = b 

pq 

It  can  be  easily  verified  t az,  .sing  (3.15) 
and  (3.4),  all  coefficients  can  be  computed.  By  means  of 
a simple  computer  program,  therefore,  we  can  compute  the 
exact  weight  that  any  measurement  has  in  the  estimation  of 
X. 

Let's  summarize  the  steps  to  compute  the  error 
variance  at  (m,n) : 

(1)  To  find  each  coefficient  b use  equation  (3.4) 

p >q 

A 

with  the  values  of  y given  by  (3.15).  The  value  of  X(m,n) 
will  be  b . The  initial  conditions  are: 

p*q 

A 

X(i,j)  = 0 for  i = (p-1)  or  j = (q-1) 


. 2 
(2)  Use  equations  (3.11)  and  (3.12)  to  compute  a . 
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(3)  Use  equations  (3.13)  and  (3.14)  to  compute  b. 

(4)  Compute  error  variance  by  (3.8). 

C.  ANOTHER  RECURSIVE  FILTER 

In  this  section  we  will  introduce  another  sub-optimum 
recursive  filter  having  the  same  structure  as  those  in  [3-5] 
It  is  similar  to  [5] , but  it  computes  the  gains  without 
approximations . 

The  objective  is  to  find  a way  to  calculate  the  gains, 
such  that  the  variance  of  the  estimation  error  is  minimized. 
Let's  call  the  covariance,  between  X and  its  prediction 

/v  2 

X , b ; and  the  variance  of  tha  prediction  a 
P P P 

Using  equation  (3.6): 


b = E (XX)  = (l-K)E(XX  ) + KE(Xy) 

b'' 


but: 


bp  - E(XXp) 


E (Xy)  = E(X2)  « a2 


Thus : 


b = (l-K)b  + Ka‘ 
P 


(3.16) 


Squaring  equation  (3.6)  and  taking  the  expected 


value : 


E(X  ) = ( 1-K)  E (Xp  ) + K E(y‘)  + 2K ( 1-K) E ( yXp) 

but: 

2 ~ 2 
°P  = E(Xp2) 

E (y2)  = a2  + r 

E(yXp)  = E(XXp)  = bp 

Therefore: 

a 2 = ( 1-K) 2a  2+K2 (c2+r)+2K(l-K)b  (3.17) 

e p p 

Substituting  equations  (3.16)  and  (3.17)  in  (3.8): 

P (m, n)  = (1-K) 2 (o2+a  2-2b  ) + K2r  (3.18) 

P P 

2 

Observe  in  (3.18)  that  c and  b are  independent  of 

P P 

the  gain  K(m,n),  but  they  are  functions  of  "past"  gains. 

2 

Therefore  we  can  find  K(m,n) , as  a function  of  a , b , 

P P 

2 

a and  r,  such  that  the  error  variance  is  minimized. 

Differentiating  (3.18)  with  respect  to  K and  equating 
to  zero: 

2 2 

0 + °n  " 2br. 

K (m,  n)  = 2 (3.19) 

eT  + a ^ - 2b  + r 


The  variance  of  the  prediction  error  is : 


P (m,n)  = E [ (X-Xp) 2 ] 


a2  + a 2 - 2b  (3.20) 

P P 


Substituting  in  (3.19): 


K (m,n) 


(3.21) 


Equation  (3.20)  is  intuitively  appealing,  because  it  is 
the  same  as  that  of  the  one  dimension  scalar  Kalman  filter. 
Substituting  (3.19)  into  (3.18): 


PMIN 


= Kr 


(3.22) 


Therefore,  with  the  gain  calculated  by  (3.21),  the 

minimum  error  variance  is  given  by  (3.22). 

Observe  that  P (m,n)  is  independent  of  K(m,n) , and, 
P 

therefore,  can  be  computed  using  "past"  gains.  Its 
calculation  is  quite  similar  to  that  of  P(m,n). 


X (m,n)  = ATY 

P P 


(3.23) 


where : 


T 

A = [au  a12  ...  apl  ...  a^^J 
Y T = [y (1, 1)  y (1, 2)  ...  y(p,l)  ...  y(m,n-l)]. 
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The  variance  of  X is: 

P 


°p  = 


(3.24) 


where  is  the  autocorrelation  matrix  of  the  measurements 
P 

-P‘ 

The  covariance  b (m,n) , similarly  to  b(m,n),  is  given 

tr 


bP  - * "xy 


(3.25) 


where  the  vector  R^Y  is  the  correlation  between  X and  Y 

P E 

The  coefficients  A are  calculated  identically  to  B, 

but  using  the  predictor  equation  (3.7)  instead  of  (3.4). 

A summary  of  gain  calculation  follows: 

(a)  Use  equations  (3.15)  and  (3.7)  to  find  the 

coefficients  A.  The  initial  conditions  for  (3.7)  are: 


X(i,j)  = 0 for  i = (p-1)  and  j = (q-1) . 

2 

(b)  Use  equation  (3.24)  to  find  a 

P 

(c)  Use  equation  (3.25)  to  find  b . 

(d)  Use  equation  (3.20)  to  find  P , and  compute 
K(m,n)  by  (3.21) . 

D.  NON- RECURSIVE  FILTER 

The  optimum  recursive  filter  must  have  the  same 
performance  of  the  non-recursive  Wiener  filter,  provided 
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the  data  set  is  the  same.  Since  two-dimensional  recursive 
filters,  like  those  in  [3-5],  are  not  optimum,  we  want  to 
find  how  good  they  are. 


The  minimum-mean-square-error  (MMSE)  in  the  estimation 

of  X(m,n),  using  the  measurements  in  region  Y (see 

m,  n 


figure  3.2),  is  computed  as  follows. 
Define  the  vector 


E ( YY 1 ) = RyY 


E (XY)  = R 


P(m,n)  = a2  - 2HTRXY  + f^R^H 


(3.26) 


Differentiate  (3.26)  and  equate  to  zero: 


H = RyY  Rj 


(3.27) 


Substitute  (3.27)  into  (3.26): 


E.  PERFORMANCE  OF  THE  FILTERS 

In  this  section  we  present  the  results  of  the  comparison 
between  filters  [3],  [5]  and  the  one  introduced  in  Section  C, 
as  well  as  the  optimum  interpolator  of  Section  D. 

First  of  all,  we  have  compared  filters  [5]  and  that  of 
Section  C,  since  they  are  essentially  the  same,  but  the 
former  uses  an  approximation.  This  comparison  was  accom- 
plished by  computing  the  gain  and  error  variance  of  both 
filters  for  several  situations.  The  result  was  that 
they  presented  practically  the  same  values;  therefore 
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we  concluded  that  the  approximation  made  in  [5]  is  very 
good.  This  conclusion  is  important,  because  it  allows  the 
use  of  simple  recursive  equations,  developed  in  [5] , for 
gain  computation. 

The  comparison  of  filters  [3]  and  [5]  against  the 
optimum  estimator  is  presented  in  table  3.1  and  figure  3.3. 

In  these  results  the  coefficient  of  correlation  is  0.9, 
in  both  directions,  and  the  variance  of  X is  1.0.  The 
variance  of  noise  is  in  the  range  0.01-1.0,  therefore  the 
signal- to-noise  ratio  goes  from  0-20  dB. 

From  figure  3.3  and  table  3.1,  we  see  that  filter  [3] 
is  better  than  [5] , but  both  are  quite  close  to  the  optimum 
filter.  In  the  worst  case,  the  error  variance  of  [3]  is 
6.5%  greater  than  the  optimum,  and  [5]  is  15%.  Also,  the 
error  variance  of  [5]  is  8%  greater  than  that  of  filter 
[3].  It  was  observed,  from  other  results,  that  [3]  and 
[5]  approach  the  optimum  for  lower  coefficients  of  corre- 
lation and/or  low  noise. 

Another  experiment  was  performed  to  change  the  gain  to  see 
what  MMSE  can  be  achieved  with  such  a filter  structure. 

Figure  3.4  and  table  3.2  show  that  filter  [3]  is  quite  close 
to  the  minimum  error  variance,  although  it  is  not  exactly 
the  minimum  as  can  be  verified  in  figure  3.5  and  table  3.3. 
The  conclusion  is  that  filter  [3]  is  better  than  [5] , but 
it  is  not  the  best  that  can  be  done.  As  a matter  of  fact 
the  present  method  of  comparison  can  always  be  used  to 
compute  the  best  gain,  although  it  is  computationally  complex 
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The  reason  why  filter  [3]  is  better  than  [5]  is  that, 

although  it  assumes  ( incorrectly)  that  the  estimation  error 

is  uncorrelated  with  the  measurements  in  region  Y . , 

m-l,n-l 

the  gain  is  chosen  such  that  the  error  is  uncorrelated  with 

y(m,n).  On  the  other  hand,  the  error  of  filter  [5]  is  not 

\ 

uncorrelated  with  y(m,n),  although  it  minimizes  the 
up-dated  error  variance  at  each  point.  What  happens  is  that 
the  value  of  such  minimum  is  not  only  dependent  on  the 
previous  error  variances,  but  also  on  the  covariances  of 
the  "past"  errors,  since  they  are  correlated. 

The  important  conclusion  is  that  the 
recursive  filters  [3],  [5]  and  the  one  of  section  C,  although 
not  optimum,  are  quite  close  to  optimality.  Filter  [3]  is 
the  best  and  simple  enough  to  be  extended  to  three  dimensions, 

0 

in  order  to  exploit  recursively  the  correlation  in  time. 

This  will  be  done  in  the  next  chapter. 
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Table  3.1 


Noise 

Variance 

Optimum 

.... 

Filter  [ 3 J 

Filter  [5] 

0.01 

0.0083729 

0.0084204 

0.0083936 

0.04 

0.0262313 

0.0269059 

0.0266656 

0.09 

0.0477314 

0.0497128 

0.0494879 

0.16 

0.0706634 

0.0743120 

i 

0.0747814 

0.25 

0.0940204 

0.0994702 

i 

0.1014409 

0.36 

0.1172805 

0.1245305 

0.1287582 

0.49 

0.1401542 

0.1491169 

0.1562232 

0.64 

0.1624769 

0.1730206 

0.1834437 

0.81 

0.1841581 

0.1961401 

0.2101098 

1.00 

0.2051533 

0.2184254 

0.2359806 

Picture  Variance  = 1.0 
Picture  Correlation  = 0.9 


I 

I 

■ 

! 
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Gain 


Error  Variance 


0.07982749 

0.26742077 

0.09579301 

0.24629241 

0.11175853 

0.23272502 

0.12772405 

0.22441578 

0.14368957 

0.21995461 

0.15965503 

0.21842539 

Filter  [3] 

0.17562050 

0.21918494 

0.19158608 

0.22175616 

0.20755148 

0.22577757 

1 

0.22351706 

0.23097461 

0.23948258 

0.23713821 

Filter  [5] 

Picture  Variance  = 1.0 
Noise  Variance  = 1.0 
Picture  Correlation  = 0.9 


ERROR  VARIANCE 


0.29999995 

0.30999994 

0.31999993 

0.32999992 

0.33999991 

0.34999990 

0.35999990 

0.36999995 

0.37999994 

0.38999993 


0.09970987 

0.09923345 

0.09892100 

0.09875959 

0.09873748 

0.09884429 

0.09907085 

0.09940886 

0.09985095 

0.10039049 


Filter  [33 


Best  Gain 


0.39999992 


0.10102153 


Filter  [5] 


ERROR  VARIANCE 


*r  * 


IV.  THREE-DIMENSIONAL  RECURSIVE  FILTER 


A.  INTRODUCTION 

The  two-dimensional  recursive  filters,  presented  in 
Chapter  III,  use  only  the  spatial  correlation  between 
pixels  in  order  to  estimate  the  gray  levels  of  noisy  pic- 
tures. In  this  chapter  we  introduce  a recursive  filter 
that  takes  advantage  of  both  correlations  in  space  and 
in  time,  when  one  has  a group  of  pictures  sequenced 
in  time.  Experimental  evidence,  presented  in  [1] , shows 
that  television  pictures  are  correlated  in  time.  This  is 
also  an  intuitive  result,  since  knowing  one  frame  we  often 
can  guess  the  next. 

The  filter  developed  here  is  an  extension  of  the  two- 
dimensional  recursive  filter  [3],  since  it  is  the  best, 
according  to  the  results  of  Chapter  III.  In  any  case,  its 
performance  will  be  analyzed,  in  order  to  evalute  the 
improvement  resultant  in  using  the  time  correlation. 


B.  FILTER  DESIGN 

The  filter  will  be  developed  under  the  following 
conditions : 

(a)  The  time-frame  is  modeled  as  an  homogeneous  random 
field  with  zero-mean  (or  known  mean)  and  autocorrelation 
function  R(i,j,k),  where  i,j  and  k are  the  distances  between 
pixels  in  the  vertical,  horizontal  and  time  coordinates, 
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X.  * 


respectively.  The  autocorrelation  is  like  the  one  intro- 
duced in  Chapter  II,  with  kernels  given  by  equation 
(2.26) . 


(b)  The  dynamic  model  is  given  by  the  partial  difference 
equation: 


X(m,n,t)  = A X + W(m,n,t) 


(4.1) 


where  A is  a row  vector  of  coefficients  and  X is  a column 
vector  of  adjacents.  The  modeling  error  W(m,n,t)  is  white 
noise  and  uncorrelated  with  the  X's  in  region  ^mnt  (see 
figure  4.1). 

In  order  to  make  clear  the  present  derivation,  let's 
assume  a specific  form  for  the  autocorrelation,  say,  first 
order  in  all  dimensions.  The  development  remains  valid  for 
the  k-order  kernels  given  by  equation  (2.26). 

The  autocorrelation  is: 


R(i, j ,k)  = a p 


n N 


i,j,k  — 0,_ 1,-2, 


(4.2) 


The  coefficients  are: 


A = [pv  ph  pfc  “PvPh  “,pvpt  iphpt  pvphpt^  (4.3) 


The  adjacents  (see  figure  4.2)  are: 


[X(m-l,n,t)  X(m,n-l,t)  X(m,n,t-1)  X(m-l,n-l,t) 


X(m-l,n, t-1)  X(m,n-l,t-l)  X(m-1 , n-1 , t-1) ] 


The  variance  of  W is 


(4.4) 


E[W2 (m,n,t) ] = Q = a2  (l-pv2)  d-Ph2) (1-Pt2) 


(c)  The  pictures  are  contaminated  with  zero-mean  white 
noise,  uncorrelated  with  the  X's. 


y(m,n,t)  = Xtm/n/t)  + v(m,n,t) 


The  noise  variance  is: 


(4.6) 


E(v  (m,n, t) ) = r 


1 . Filter  Structure 

Since  we  want  to  generate  a random  field  X as  close 
as  possible  to  X,  let's  choose  a structure  similar  to  the 
dynamic  model  (4.1): 


X(m,n,t)  = B X + K (m,n, t) y (m,n, t) 


(4.7) 


In  what  follows  we  will  drop  the  argument  (m,n,t) , 
where  it  could  not  cause  confusion. 
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Observe  that  X can  be  written  as  a linear  combina- 
tion of  all  measurements  within  region  Y defined  in 

mnt 

figure  4.3.  The  optimum  recursive  filter  must  have  the 
gains  such  that  the  weights  of  the  measurements  are  exactly 
the  same  as  for  the  non-recursive  estimator,  constrained 
to  the  same  data  set  Y The  criteria  of  optimality  is 

minimization  of  the  Mean  Square  Error  (MSE) . It  will  be 
seen  that  it  is  impossible  to  have  an  optimum  recursive 
filter,  constrained  to  the  dynamic  model  (4.1). 

2 . Orthogonality  Principle 

A necessary  condition  of  optimality  is  that  the 
estimation  error  be  uncorrelated  with  the  data  set.  Let's 
apply  this  condition  to  find  the  coefficients  B and  K. 

The  estimation  error  is : 


The  orthogonality  condition  is : 


for  all  (p,q,r)  in  ^mnt  (see  figure  4.3) 

First,  let's  apply  (4.9)  to  the  measurement  y(m,n,t) 


/v  /v 

E (ey ) = E[(X-x)y]  = E [y (X-BX-Ky ) ] 


= E [Xy-Ky2-ByX] 


Figure  4 . 3 

Definition  of  the  region  Ym^t 


X 


With  such  a choice,  the  second  term  of  (4.12) 


<%  m 


becomes : 


(l-K)E[W(m,n,t)y (p,q,r) ] 


t 


Where  we  have  used  (4.1),  and  this  term  vanishes 
because  W(mrn,t)  is  uncorrelated  with  all  measurements  in 
region  X Substituting  (4.13)  into  (4.10): 


+ b_ 


K(m,n,t)  = 


+ b + r 
P 


(4.14) 


where  we  have  defined: 


X (m,n, t)  = AX 

(4.15) 

/N 

b (m,  n , t)  = E (XX  ) 


Equation  (4.15)  seems  to  be  an  adequate  choice  for 
the  first-step  predictor,  also  (4.14)  is  an  intuitive  result 
for  the  gain,  since  it  is  unity  for  zero  noise  (r=0)  and 
decreases  when  noise  increases. 

Unfortunately,  the  first  term  of  (4.12)  remains  and 
we  could  not  force  it  to  vanish.  Choosing  K by  (4.14)  we 
have  forced  the  estimation  error  to  be  uncorrelated  with 
y(m,n,t),  but  it  is  not  correct  to  conclude,  by  induction, 
that  the  estimation  errors  of  the  adjacents  (vector  e)  are 
uncorrelated  with  the  other  measurements. 


^ 


iff 

I 
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Take,  for  example,  the  error: 


e(m,n,t-l)  = X(m,n,t-1)  - X(m,n,t-1) 

The  estimate  didn't  use  any  measurements  of  frame 
t,  therefore  the  error  can  not  be  uncorrelated  with  them. 

As  a matter  of  fact,  the  errors  might  only  be 
uncorrelated  with  the  measurements  within  region  Y . . ^ . . 

In  reference  [3],  for  the  two-dimensional  case,  it  is 
stated  that  the  error  is  uncorrelated  with  all  measurements 
in  a similar  region.  Unfortunately,  that  is  not  correct. 

The  estimation  error,  in  the  two  or  three  dimensional  cases, 
is  uncorrelated  only  with  y(m,n,t),  because  K was  computed 
under  such  conditions,  and  with  the  measurement  y( 1,1,1). 

This  can  be  seen  from  (4.12)  and  (4.13): 


E[ey(p,q,r)]  = (1-K)E[A  e y(p,q,r)j 


for  (p, q , r ) f (m,n,t) 


(4.16) 


Initializing  the  filter  with  the  mean  value  of  X 
(=0),  it  can  be  verified,  using  (4.16),  that  e(m,n,t)  is 
uncorrelated  with  y( 1,1,1). 

Although  equation  (4.14)  is  not  the  optimum  choice 
» for  the  gain,  we  have  shown  in  Chapter  III  that  it  was  the 

\ best  in  the  two-dimensional  case,  therefore  we  hope  it  is 

also  good  in  three  dimensions.  This  will  be  verified  at 
the  end  of  this  chapter. 


Using  (4.7)  and  (4.13),  the  filter  equation  becomes: 


X (m, n, t)  = (l-K)A  X + Ky 


(4.17) 


Alternatively,  using  (4.15): 


X (m,n, t)  = (l-K)Xp  + 


(4.18) 


3.  Gain  Computation 

V 

Let's  develop  a recursive  equation  to  compute  b 

h 

and,  therefore,  the  gain.  Define  a function: 


F (m,n, t)  = E[X(m,n,t)X(p,q,r) ] 
pqr 


Using  (4.17) : 


F = [ 1-K (p , q, r ) ] A F + K(p,q,r) R (m-p , n-q , t-r ) (4.19) 

pqr  


where : 


T 

— = ^Fp-l,q,r  Fp,q-l,r  Fp,q,r-1  Fp-l,q-l,r  p-l,q,r-l 


Fp,q-l,r-l  Fp-l,q-l,r-l^ 


(4.20) 


Initializing  with  X equal  to  the  mean  of  X,  outside 
the  picture,  the  initial  conditions  for  F are: 


96 


F (m,n,t)  = 0 for  p,q  or  r = 0 

P*4»r 

Now,  using  (4.15)  and  (4.19): 

b_(m,n,t)  = A F (4.21 

P p 

where  is  the  value  of  F for  (p,q,r)  = (m,n,t) . 

Summary  of  gain  calculation: 

(a)  Use  (4.19)  to  compute  each  component  of  the  vector 

~P* 

(b)  Use  (4.21)  to  compute  b (m,n,t). 

b' 

(c)  Compute  K(m,n,t)  by  (4.14) 

Observe  that  the  calculation  of  each  gain  K(m,n,t) 
requires  scanning  the  time-frame  from  (1,1,1)  through  the 
pixel  (m,n,t).  Fortunately,  as  will  be  seen,  after  a few 
frames  the  gain  reaches  steady  state,  therefore  the  compu- 
tation can  be  reduced  to  the  top  left  corner  of  a few 
frames . 

C.  PERFORMANCE  OF  THE  FILTER 

The  method  used  in  Chapter  III,  for  the  two  dimensional 
case,  can  be  easily  extended  to  three  dimensions.  The 
variance  of  the  estimation  error  is : 

= E[ (X  - X) 2] 

= E (X2)  + E (X2 ) - 2E (XX) 


P (m,n,t) 


or: 


P(m,n,t)  = a2  + a 2 - 2b 

e 


(4.22) 


where: 


2 2 
E (X  j = a 


"2  2 
E(XZ)  = 


E (XX)  = b 


The  covariance  b can  be  calculated  using  b and  K: 

P 


b (m,n, t)  = E (XX)  = (l-K)E(XXp)  + KE(Xy) 


b(m,n,t)  = (l-K)b  + Ka‘ 


(4.23) 


To  compute  the  variance  of  the  estimator,  write  X as 
a linear  combination  of  the  measurements  in  region  Y_  . : 


X(m,n,t) 


(4.24) 


where  c is  a column  vector  of  coefficients , function  of 
the  gains,  and  Y is  a column  vector  containing  all  measure- 
ments in  Y 

mnt 

A 

The  variance  of  X is : 
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2 T 

ae  = C V 


(4.25) 


where  Ry  is  the  correlation  matrix  of  the  vector  Y: 


Ry  = E(YY) 


The  key,  in  this  computation,  is  the  method  of  finding 
the  coefficients  c.  This  can  be  easily  accomplished 
using  the  filter  equation  (4.17)  with: 


y(i,j,k)  = 


1 (i,j,k)  = (p,q,r) 


(i,j,k)  j-  (p,q,r) 


X(i,j,k)  = 0 for  k = (r-1)  or 

i = (p-1)  or 
j = (q-1) 


Initializing  (4.17)  as  above,  and  starting  at  (p,q,r) , 
the  resultant  value  at  (m,n,t)  will  be  the  coefficient  of 
y(p,q,r)  in  the  estimate  X(m,n,t).  In  this  way  all  coeffi- 
cients can  be  found  and,  therefore,  the  error  variance. 

We  are  also  interested  in  the  first-step  predictor 
given  by  (4.15).  Its  error  variance  is: 


P (m,n,t)  = E[(X  - Xpr] 


(4.26) 
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It  can  be  related  with  P(m,n,t)  using  (4.18)  and  (4.26) 


X - X = X - (1-K) X - Ky 

P 


= (1-K) (X  - X ) - KV 

P 


Squaring  and  taking  expectations : 


P = (1-K)2P  + K2r 

P 


Pp(m,n,t) 


(P  - K r) 


(1-K) 


(4.27) 


D.  RESULTS 


In  figures  (4.4)  through  (4.9)  some  results  of  gain 
calculations  are  shown.  For  the  first  frame  the  gain  is 
the  same  as  in  [3]  for  two  dimensions.  From  these  results 
we  can  observe  that  the  gain  reaches  a steady-state  value 
very  fast,  at  about  frame  number  4.  Also,  for  the  same 
frame,  the  steady-state  is  reached  at  about  pixel  (4,4). 

That  is  an  interesting  result,  since  it  simplifies  substan- 
tially gain  calculations.  Also  we  can  think  of  using  a 
constant  gain  for  the  whole  time-frame,  or  a few  gains 
for  the  first  4 frames. 

In  tables  4.1  through  4.3  a comparison  is  shown 
between  the  three-dimensional  filter,  estimate  and  prediction, 
and  the  two-dimensional  recursive  filter.  It  can  be  seen 
in  table  4.1  that  exploitation  of  time  correlation  can  be 
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quite  advantageous  in  cases  where  time  correlation  is 
very  high,  compared  with  spatial  correlations.  In  this 
example,  the  estimation  error  variance  is  reduced  by  about 
33%,  due  to  the  use  of  time  correlation.  In  table  4.2 
a comparison  with  equal  amounts  of  correlation  in  time  and 
space  is  shown.  In  this  example  the  improvement  was  about 
10%.  It  can  be  seen  in  table  4.3  that  exploitation  of  time 
correlation  is  not  so  advantageous  when  the  images  are 
more  correlated  in  space  than  in  time.  In  this  example 
the  improvement  was  only  around  6%. 
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SIGNAL/NOISE  = 10  dB 

FRAMES  CORRELATION:  p = p,  = = 0.8 

V n t 


i e 


Figure  4.5 


K 


Gains  for  the  Time-Frame 
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SIGNAL/NOISE  = 30  dB 

FRAME  CORRELATION:  pv  = ph  = 0 . 9 


Figure  4 . 6 

Gains  for  the  First  Frame 


Table  4.1 


Noise 


Error 

Variance 

Three-Dimensions 

Two-Dimensions 

Estimation 

Prediction 

Estimation 

Prediction 

.236 

.293 

.357 

.554 

.210 

.267 

.318 

.527 

.178 

.236 

.271 

.493 

.140 

.197 

.212 

.450 

.089 

.144 

.132 

.387 

Picture  Variance  = 1.00 
Spatial  Correlation  = 0.70 
Time  Correlation  = 0.95 


Table  4.2 


Noise 

Variance 


Error 

Variance 

Three-Dimensions 

1 

Two-Dimensions 

Estimation 

Prediction 

Estimation 

Prediction 

.269 

. 351 

.301 

.427 

.240 

.324 

.268 

.400 

.206 

.293 

.229 

.367 

.162 

.251 

.180 

.324 

.103 

.191 

.114 

.263 

Picture  Variance  = 1.0 
Spatial  Correlation  = 0.8 
Time  Correlation  = 0.8 


Error  Variance 


Three-Dimensions 

Two- Dimens ions 

Estimation 

Prediction 

Estimation 

Prediction 

.280 

.376 

.301 

.427 

.250 

-349 

.268 

.400 

.214 

.317 

.229 

.367 

.168 

. 274 

.180 

.324 

.107 

.213 

.114 

.263 

- J 

Noise 

Variance 


1.0 

0.8 

0.6 

0.4 

0.2 


Picture  Variance  = 1.0 
Spatial  Correlation  = 0.8 


Time  Correlation  = 0.7 


A.  INTRODUCTION 

In  this  chapter  a new  class  of  image  filters  is  intro- 
duced, called  hybrid  filters.  The  recursive  filters, 
presented  in  Chapters  III  and  IV,  utilize  only  part  of  the 
data  set,  arbitrarily  defined  as  "past"  measurements.  Since 
the  whole  picture  is  often  available,  an  optimal  smoother 
might  be  sought  in  order  to  utilize  all  the  information 
available.  For  the  two-dimensional  recursive  filters  of 


Chapter  III,  the  smoother  would  be  the  optimum  combination 
of  the  four  estimates  obtained  by  scanning  the  picture 
starting  at  each  corner.  The  first  difficulty  is  that 
those  filters  are  not  optimum  resulting  in  a sub-optimum 
smoother.  Second,  the  smoother  would  require  scanning  the 
picture  four  times,  therefore  increasing  substantially  its 
complexity  of  implementation. 

The  hybrid  filter  introduced  here  is  a smoother  that 
combines  optimally  the  estimate  of  the  recursive  filter 
(two  or  three  dimensions)  with  an  arbitrary  set  of  "future" 
measurements.  It  will  be  applied  in  picture  enhancement 
and  compared  against  the  recursive  filter  as  well  as  the 
non-recursive  filter  presented  in  [6] . Some  examples  of 
hybrid  filters  designed  to  predict  the  pixel  gray  level  will 
also  be  presented.  These  filters  will  be  applied  in 
Chapter  VI  for  purposes  of  target  detection. 
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The  hybrid  filter  is  particularly  useful  in  some 
applications  where  the  pictures  are  not  from  a fixed 
scenario.  For  example,  the  camera  has  a movement  in  order 
to  follow  some  target.  In  such  cases,  the  background  may 
change  so  much  that  past  frames  don't  carry  enough  informa- 
tion about  the  present  frame.  In  this  case,  the  three- 
dimensional  recursive  filter  can't  be  applied,  but  a hybrid 
filter  can  be  designed  using  recursively  the  measurements 
in  the  present  frame,  and,  non-recursively , the  "crude" 
measurements  in  the  previous  frame. 

B.  TWO-DIMENSIONAL  HYBRID  FILTERS 

In  this  section  we  present  two  examples  of  hybrid 
filters  in  two-dimensions.  One  is  a smoother  designed  to 
enhance  pictures  contaminated  with  additive  white  noise. 

The  other  is  a filter  designed  to  predict  the  pixel  gray 
level,  based  on  noisy  measurements  of  previously  scanned 
pixels.  Here  we  assume  the  picture  is  scanned  row  by  row 
from  top  to  bottom. 

1.  Two-Dimensional  Hybrid  Smoother 

The  recursive  filters  presented  in  Chapter  III 
utilize  the  measurements  in  region  Y (see  figure  3.2) 

A 

to  have  an  estimate  X(m,n).  To  improve  this  estimate  we 
are  going  to  use  also  a selected  set  of  measurements  in 
the  neighborhood  of  (m,n) , say  Y,  and  combine  optimally 

A A 

with  X(m,n),  in  order  to  have  a better  estimate  X^(m,n). 
Let's  choose  for  Y the  set  of  five  neighbors  shown  in 
figure  5.1. 


Ill 


A 


P(m,n)  = E [ (X  - hT2.)  2 ] (5.2) 

Differentiating  (5.2)  with  respect  to  h and  equating 
to  zero 


where 


Rgg  - E(aaT) 

KXg  = E(xa) 

The  measurement  equation  is 


y(m,n)  = X(m,n)  + v(m,n)  (5.4) 


where  v is  zero-mean  white  noise,  uncorrelated  with  X, 
having  variance 

2 

E [v  (m,n) ] = r 


The  gray  level  X is  a zero-mean  (or  known  mean) 
homogeneous  random  field  with  autocorrelation 

R ( i , j ) = E [ X (k , £) X (k+i , £+ j ) ] (5.5) 

The  elements  of  the  autocorrelation  matrix  R , that 

gg 

A 

don't  include  X(m,n),  are  easily  found  using  equations 

A 

(5.4)  and  (5.5).  To  find  the  elements  which  include  x, 
observe  that  in  Chapter  III  a method  was  developed  to  cal- 
culate error  variance  of  recursive  filters  [3,5],  where 
E(X^)  is  also  computed  as  well  as  E(XX).  It  can  be 
verified  that  the  same  method  can  be  used  for  higher  order 
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V 


recursive  filters,  though  in  this  derivation  we  restrict 
ourselves  to  recursive  filter  [3]  . 


i 


t 

m 

o 


i 


Therefore,  using  (3.11)  and  (3.13): 

2 

E(X^)  = 

(5.6) 

A 

E (XX)  = b 


It  can  be  verified  that  the  covariances  between 

A 

X and  the  measurements  Y can  be  expressed  as  a function 
of  b given  by 


E [Xy (m+l,n-l) ] = PvPhb 

/v 

E [Xy (m+1 , n) ] = pyb 

A 

E[Xy (m+l,n+l) ] = p^p^b  (5.7) 

/\ 

E [Xy  (m,n+l)  ] = phb 

A 

E [Xy  (m-l,n-l)  ] = PvPhb 

Therefore  we  have 
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a.  Performance  and  Comparison 


The  error  variance  of  the  hybrid  smoother  is 
calculated  substituting  (5.3)  into  (5.2): 


P (m, n) 


(5.9) 


At  this  point  it  is  interesting  to  compare  the 
hybrid  smoother  and  the  recursive  filter  [3] , since  we 
know  how  to  calculate  the  error  variance  for  both.  It  is 
also  interesting  to  include  in  such  comparison  the  non- 
recursive filter  presented  in  [6] . This  filter  uses  a 
small  window  of  3x3,  shown  in  figure  5.2,  to  estimate  the 
pixel  gray  level  at  the  center,  minimizing  the  mean-square- 
error.  The  hybrid  smoother  must  have  superior  performance, 
since  it  utilizes  the  same  measurements  of  [3]  and  [6] 
together.  However,  we  can't  say  which  is  better  whether 
[3]  or  [6],  since  they  use  distinct  data  sets.  Although 
filter  [3]  uses  a larger  data  set  than  [61 , the  former 
uses  all  the  closest  adjacents  to  the  estimated  pixel 
(m,n)  . 

In  table  5.1  some  numerical  results  are  shown 

to  help  the  comparison.  As  expected,  the  hybrid  filter 

presented  the  best  performance  with  error  variance  about 

25%  smaller  (for  p = 0.94)  than  filter  [6].  The  non-recursive 

filter  [6]  presented  better  performance  than  the  recursive 

filter  [3] . The  error  variance  of  the  optimum  interpolator 

constrained  to  the  observations  in  Y is  also  included  in 

m , n 
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Table  5.1 


Error  Variance 

Correlation 

Hybrid 

Filter  [6] 

Filter  [3] 

Optimum 

Interpolator 

0.60 

0.334 

0.346 

0.400 

0.399 

0.70 

0.279 

0.298 

0.357 

0.355 

0.80 

0.217 

0.245 

0.300 

0.294 

0.90 

0.145 

0.182 

0.218 

0.205 

0.92 

0.128 

0.167 

0.197 

0.181 

0.94 

0.112 

0.152 

0.173 

0.155 

0.96 

0.111 

0.136 

0.146 

0.124 

Picture  Variance  = 

L.O 

j 

1 

i 

j Noise  Variance  = 1.0 

I 

i 

this  table.  It  is  interesting  to  see  that  the  use  of  just 
5 non-causal  observations  by  the  hybrid  filter  were  enough 
to  fully  compensate  the  suboptimality  of  the  recursive 
filter  [3] . Observe  that  the  9-point  non-recursive  filter 
[6]  also  presented  better  performance  than  the  optimum 
interpolator,  even  for  spatial  correlation  as  high  as  0.94. 
This  indicates  that  most  of  the  information  about  X(m,n) 
resides  in  its  nearest  neighbors.  This  motivates  the  deriva- 
tion of  simple  filters  like  the  hybrid  filter  introduced 
here.  It  is  much  simpler  than  the  optimum  recursive  filter 
[9]  with  insignificant  loss  in  performance.  It  also  has 


I 


the  same  order  of  complexity  as  the  non-recursive  filter 
[6],  but  superior  performance.  Complexity  here  means  the 
number  of  multiplications,  which  is  9 for  the  non-recursive 
filter  and  10  for  the  hybrid  filter. 

The  hybrid  filter  may  be  seen  as  a way  of 
artificially  increasing  the  window  of  filter  [6] . 
b.  Experimental  Results 

In  this  section  we  present  the  results  of  some 
experiments  where  the  filters  were  used  in  picture  enhance- 
ment. In  this  experiment,  the  picture  of  figure  5.3  was 
contaminated  with  white  noise.  Assuming  knowledge  of  the 
noise  variance,  the  picture  was  modeled  by  the  method 
presented  in  Chapter  II.  With  this  model  and  noise  vari- 
ance, we  designed  the  three  filters:  hybrid,  recursive 
and  non- recursive . The  theoretical  performance  was  computed 
for  each  filter,  as  in  the  last  item.  A signal-to-noise 
ratio  (S/N)  was  defined  as  the  ratio  between  the  variance 
of  the  original  picture  (figure  5.3)  and  the  variance  of 
the  noise.  A processing  gain  was  defined  as  the  ratio  of 


Table  5.2  summarizes  the  results.  It  can  be 
seen  that  the  computed  performance  is  pretty  close  to  the 
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results  of  the  experiments.  The  hybrid  is  clearly  superior 
to  the  others  and  the  non-recursive  filter  is  better  than 
the  recursive,  except  for  the  very  low  S/N  of  -3  dB,  when 
the  former  was  a little  better.  However  the  subjective 
evaluation  of  the  pictures  (figures  5.4  through  5.15)  is 
more  difficult.  In  our  evaluation  the  hybrid  filter  had 
superior  performance,  therefore  validating  the  theoretical 
results,  as  well  as  the  measurements  of  table  5.2.  However, 
the  results  of  the  non-recursive  filter,  in  our  subjective 
evaluation,  were  inferior  to  the  recursive  filter  in  all 
three  experiments,  contradicting  the  measurements  of  table 
5.2.  Probably  tne  mean-square-error  is  not  a good  figure 
of  merit  for  visual  perception,  however  there  is  no  doubt 
that  considerable  improvement  was  achieved  by  processing 
the  noisy  pictures  using  such  criteria. 

2 . Two-dimensional  Hybrid  Predictor 

Assuming  that  the  picture  is  scanned  row  by  row, 
we  may  define  "past"  measurements  as  those  pixels  already 
scanned.  Then  the  gray  level  of  pixel  (m,n)  may  be  pre- 
dicted using  "past"  measurements.  In  figure  5.16  the 
configuration  of  a hybrid  filter  designed  to  predict  X(m,n) 
is  shown.  This  filter  combines  optimally  the  prediction 

A 

Xp(m,n)  of  the  recursive  filter  and  the  measurements  Y, 
shown  in  figure  5.16.  The  result  is  an  improved  prediction 

A 

X(m,n)  that  utilizes  all  the  information  contained  in  the 

measurements  Z „ and  Y. 

mn  — 


Figure  5.10 


Image  of  figure  5.8  filtered  by  the  recursive  filter 


(SNR  = 7.8  dB) 


The  design  of  this  filter  is  similar  to  the  hybrid 
smoother  of  the  previous  item,  therefore,  we  present  below 
the  results. 

The  autocorrelation  function  used  is 


R(i,  j) 


The  predictor  is 


(5.1^ 


X(m,n)  = hT£ 


where 


r 

A 

r~ 

X (m,n) 

hr 

hi 

! y(m-l,n+l) 

1 

h2 

| y(m-l,n+2) 

h = 

h3 

y (m-2 ,n+l) 

h4  i 

y (m-2 ,n+2) 

h5  ! 

_ 

L J 

The  weights  h are  given  by 


Using  the  results  of  Chapter  III  we  can  compute 


all  the  elements  of  R _ and  R..  . Let's  call: 


1 ' m ww 


E [X  2 (m, n) ] = o 2 

P P 


(5.11) 


E[XXp]  = bp 


It  can  be  verified  that  the  cross-correlation 

/N 

between  X and  Y can  be  expressed  as  a function  of  b , 
P P 


as  follows: 


E[Xpy  (m-1 , n+1)  ] 


o o,  b 
v h p 


E[Xpy  (m-1,  n+2 ) ] 


2, 

p p,  b 

vMh  p 


(5.12) 


E [X  y (m-2  ,n+l)  ] 

P 


p 2p.b 

Hh  p 


E[X  y (m-2, n+2)  ] 

P 


2 2U 

pv  ph  bp 


Using  equations  (5.10),  (5.11)  and  (5.12),  we  can 


calculate  the  autocorrelation  matrix  R and  the  correlation 

99 


T 


■*  *■ 
o 


P p.cr 

v h 


2 2 

*Xg  Pvph  a 

P Va2 

v h 


n 2 2 2 
pv  Ph  ° 


C.  THREE-DIMENSIONAL  HYBRID  FILTERS 

In  Lhis  section  we  present  two  hybrid  filters 
applied  to  pictures  sequenced  in  time-  The  first  filter 
is  a smoother  to  be  used  in  image  enhancement  and  the 
other  is  a predictor.  The  design  of  these  filters  is 
quite  similar  to  those  of  the  previous  section. 

1.  Three-dimensional  Hybrid  Smoother 

The  filter  structure  is  shown  in  figure  5.17. 

The  recursive  filter  is  that  developed  in  Chapter  IV.  The 

A 

estimate  X(m,n,t)  of  the  recursive  filter  is  optimally 
combined  with  the  set  of  measurements  Y,  in  order  to  have 
a better  estimate  X^(m,n,t).  The  design  is  quite  the  same 
as  in  the  two-dimensional  case,  thus  we  present  only  the 
results . 

The  autocorrelation  used  for  the  time-frame  is 


R(i,j,k) 


(5.14) 
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The  estimate  is 


~ T 

X(m,n,t)  = h £ 


where : 


X(m,n,  t) 

hl 

y (m+l,n-l, t) 

h2 

y(m+l,n,t) 

h3 

i y (m+l,n+l, t) 

| 

h4  i 

y (m, n+1 , t) 

hs  ! 

y (m-l,n+l, t) 

h = 

h6  I 

y (m+1 , n-1 , t-1) 

h7  ! 

j y (m+1, n, t-1) 

h8 

y (m+1, n+1, t-1) 

h9  ; 

: 

y (m, n+1, t-1) 

hio  i 

y (m-1 , n+1 , t-1) 

1 

H 

i — 1 

-C 

1 

The  weights  h are  given  by 


Using  the  same  notation  as  in  Chapter  IV  , we  have 


E [X2 (m, n, t ) ] * a 2 

e 


(5.15) 


As  before,  the  correlation  between  X and  the  measure 
ments  Y can  be  expressed  as  a function  of  b as  follows 

E [X(m,n, t)y (m+i,n+j ,t+k) ] = Pv^  PjJ  ^ b (5.16) 


where 


(i,j,k)  = (1,-1, 0) , (1,0,0) , (1,1,0) , (0,lk0) , 

(-1,1,0)  (1,-1,-1)  (1,0,-1) , (1, 1,-1)  , 
(0,1,-1) , (-1,1,-1) 


. equations (_5._M),  (5.15)  and  (5.16),  we  can 

calculate  the  autocorrelation  matrix  R and  the  correlation 

gg 

vector  R , then  the  weights  h. 

xg  — 

2.  Three-dimensional  Hybrid  Predictor 

This  filter  combines  optimally  the  first-step 

A 

prediction  X of  the  two-dimensional  recursive  filter  [3] 
and  the  set  of  measurements  Y shown  in  figure  5.18.  Since 
its  design  is  quite  similar  to  the  two-dimensional  case, 
we  present  only  the  results. 

The  autocorrelation  function  is  given  in  (5.14) 
and  the  hybrid  prediction  is 

A _ 

X(m,n,t)  = h £ 


where : 


Xp (m,n) 

hl 

y (m-l,n-l, t-1) 

h2 

1 

! y(m-l,n,t-l) 

h3 

y (m-l,n+l, t-1) 

h4 

y (m,n-l, t-1) 

h = 

h5 

y (m,n, t-1) 

1 

h6 

y (m+l,n+l, t-1) 

h7 

! y (m+l,n-l, t-1) 

00 

i y(m+l,n,t-l) 

! hq  i 

y i 

y (m+l,n+l, t-1) 

[-10 

The  weights  h are  given  by 


Similarly  to  the  two-dimensional  case: 


(5.17) 


E[XXp]  = bp 


And  the  correlation  between  Xp 


and  Y is  given  by 


E (Xp  (m, n, t) y (m+i , n+ j , t— 1 ) ] 


p . b 

t p 


(5.18) 


where 


The  autocorrelation  matrix  Rgg  an<^  the  correlation 

vector  Rv  can  be  calculated  using  equations  (5.14) , (5.17) 

Xg 

and  (5.18),  then  the  weights  h are  computed. 

D.  COMPARISON  OF  THE  PREDICTORS 

The  predictors  developed  here  will  be  used  in  the  target 
detection  problem  addressed  in  Chapter  VI.  In  this  section 
we  present  some  numerical  results  in  order  to  evaluate  its 
performance. 

In  Tables  5.3  through  5.5  the  error  variance  of  the 
recursive  and  hybrid  predictors  for  several  situations  is 
shown.  From  these  results  we  can  see  that  the  exploitation 
of  the  time  correlation  can  be  quite  advantageous , mainly 
for  the  case  of  high  correlation  in  time  (table  5.3),  where 
the  error  variance  of  the  three-dimensional  filters  is  around 
-~rfta  n-F  rr i Tran  by  the  two~dimens ional  filters.  We  also  can 


observe  in  the  three  tables  that  Htp  -£i_J  i s 

better  them  the  3D-recursive  filter  for  observation  noise 
variances  above  0.2,  although  the  former  exploit  only  the 
time  correlation  with  the  previous  frame.  This  result  rein- 
forces the  basic  idea  of  the  hybrid  filters,  that  is  the  use 
of  all  the  neighbors  of  uhe  pixel  (m,n) , including  those 


that  are  not  causal. 


Table  5.5 


Noise 

Variance 

Error 

Variance 

three-dimensional 

two-dimens ional 

recursive 

hybrid 

recursive 

hybrid 

1.0 

.376 

.338 

.427 

.364 

00 

• 

o 

.349 

.317 

.400 

.340 

0.6 

.317 

.292 

.367 

.312 

0.4 

.274 

.260 

.324 

.277 

0.2 

.213 

.216 

.263 

.227 

Picture  variance  = 1.0 


pv  = ph  = °.8 

Pt  - 0.7 


A.  INTRODUCTION 


In  Chapter  II  we  have  developed  models  for  images  and 
for  a sequence  of  images  in  time.  Using  such  models  we 
have  constructed  two  and  three  dimensional  recursive  and 
hybrid  filters  to  smooth  out  noise  in  images.  In  this 
chapter  we  are  interested  in  the  detection  of  targets  using 
low  contrast  images.  This  problem  is  usually  referred  to 
as  background  suppression,  where  background  means  a kind 
of  texture  that  predominates  in  the  picture.  In  this  con- 
text, the  word  target  means  the  information  that  we  are 
interested  in  extracting  from  the  pictures.  It  may  be  a 
ship,  an  airplane  or  some  kind  of  texture  statistically 
different  from  that  of  the  background. 

Both  textures,  target  and  background,  are  modeled  as 
in  Chapter  II,  but  it  is  assumed  that  less  "a  priori" 
knowledge  is  known  about  the  target.  The  basic  idea  in 
detection  is  to  use  the  prediction  of  the  image  filters 
previously  developed.  Knowing  this  prediction  and  the  actual 
observation  of  the  pixel  we  develop  here,  a decision  rule 
to  decide  to  which  texture  the  pixel  belongs. 


t 

< 

■ 
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B.  TARGET  AND  BACKGROUND  REGIONS 

Figure  6.1  shows  two  distinct  regions  in  the  picture: 
background  and  target  regions.  In  the  visible  spectrum 
what  happens  is  a replacement  process,  we  have  background 
or  target.  On  the  other  hand,  in  the  infrared  spectrum 
what  seems  to  happen  is  some  combination  of  the  radiations 
from  the  target  (e.g.,  a plane)  and  the  background (e.g. , 
clouds).  To  include  these  situations,  let's  assume  that  the 
gray  level  in  the  target  region  is  given  by 


y(m,n)  = ax(m,n)  + T(m,n)  + v(m,n) 


(6.1) 


This  model  assumes  a linear  combination  of  target  and 
background,  normalized  with  respect  to  the  target,  and  also 
consider  an  additive  white  observation  noise  v(m,n), 
uncorrelated  with  both  target  and  background  textures . The 
background  texture  x(m,n)  is  modeled,  as  in  Chapter  II,  as 
a homogeneous  random  field  with  known  mean  and  autocorre- 
lation function  and  it  is  considered  independent  of  the 
target  texture.  The  target  region  is  similarly  modeled, 
but  the  amount  of  "a  priori"  knowledge,  for  most  applica- 
tions, may  be  restricted  to  the  mean  and/or  variance. 

The  background  region  is  modeled  by 


y(m,n)  = x(m,n)  + v(m,n) 


T6.2) 


U-- 


The  term  ax(m,n)  in  (6.1)  takes  into  account  a possible 
correlation  between  the  two  regions.  The  weight  a is  a 
factor  to  be  adjusted  according  to  the  specific  nature  of 
the  picture.  In  the  visible  spectrum  it  seems  reasonable 
to  assume  a replacement  process,  therefore  a = 0.  In  the 
infrared  spectrum  it  may  assume  a constant  value  at  the 
edges  of  the  target  and  zero  inside  this  region,  provided 
we  assume  that  the  sensor  has  rejected  the  background. 


C.  LIKELIHOOD  RATIO  TEST 

In  this  section  we  will  construct  a likelihood  ratio  to 
test  the  two  hypotheses  against  a threshold.  We  will  first 
consider  the  general  case  of  deciding  between  two  random 
fields,  target  and  background,  where  the  regui*^~>,'3r' priori"  ~ 
knowledge  are  the  mean  and_,aoitodorrelation  function  of  the 
background  arift-tiTe*  mean  and  variance  of  the  target  and  the 
_,_jReasurement  noise.  After,  we  will  work  some  special  cases 
and  derive  its  R.O.C.  (Receiver  Operating  Characteristic). 

The  detector  developed  here  will  be  a Newman-Pearson  detector 
[15],  where  the  threshold  is  obtained  from  the  R.O.C.  sub- 
jected to  the  constraint  that  the  probability  of  false 
alarm  is  less  than  some  desired  value. 

The  two  hypotheses  are: 

Hypothesis  H„  - The  measurement  y(m,n)  is  in  the  back- 


ground region: 


y!Ho  - yQ  = x + v 


Hypothesis  - The  measurement  y(m,n)  is  in  the  target 


region: 


Y H, 


*1 


= ax  + T + v 


The  decision  rule  will  be  applied  to  each  pixel  using 
the  observation  y(m,n)  and  the  prediction  of  the  pixel  gray 

A 

level  x (m,n) . Assume  that  the  picture  is  being  scanned 
row  by  row,  from  top  to  the  bottom.  For  the  sake  of  clarity, 
assume  that  the  past  pixels  were  correctly  identified  as 
belonging  to  the  background  region.  The  image  filters 
developed  before  (recursive  or  hybrid)  are  used  to  give  the 

A 

prediction  xp(m,n).  We  are  looking  for  a decision  rule 

A 

based  only  on  y(m,n)  and  Xp(m,n). 

Let's  define  the  likelihood  ratio: 


L(y,xr 


= In 


P(xpy |HX) 
p(xpy|H0) 


(6.3) 


where  p(xpy|’)  is  the  joint  probability  density  function, 
conditioned  to  the  hypothesis  Hq  or  . 

Applying  Bayes  rule : 


P(x  y|'Hi)  = p(x  |Hi)p(y|x  Hi)  i = 0,1 


Since  the  prediction  xp  is  independent  of  the  hypothesis: 


p'vv  = 


PUpiV 
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Lt 


J\i.” 


Substituting  these  results  into  (6.3) 


L(y , x ) = In 


p(y lxpHl) 


P(y|xpHo) 


(6.4) 


At  this  point  we  will  make  the  assumption  that  the 
conditional  probabilities  density  functions  in  (6.4)  are 
Gaussian. 

We  can  write 


y = x + v — e 
1 o p p 


y | Ht  = ax„  + v - ae„  + T 
1 1 P P 


Since  x is  known,  the  conditional  densities  are  func- 
P 

tions  of  (v-e  ) and  (v-ae  +T) , where  e is  the  prediction 

hr  hr  b' 

error.  In  figure  (6.2)  the  histogram  of  (v-e  ) is  shown  for 
the  noisy  image  in  figure  (5.8),  using  the  first-step 
predictor  of  the  two-dimensional  recursive  filter  [3] . 

This  result  shows  a curve  quite  close  to  a normal  distri- 
bution, and,  therefore,  validates  the  above  assumption  for 

A A 

p(y|x  H ) and  also  for  p(y|x  H,),  provided  the  target  is 

hr 

deterministic . 

Since  the  conditional  probabilities  density  functions 
are  assumed  Gaussian,  we  need  only  to  find  its  mean  and 
variance.  Before  proceeding,  let's  first  define  some 

A 

statistics  for  y and  x . 

P 


* 


Figure  6.2 


Histogram  of  (v-e  ) for  the  image  in  figure  (5.8) 

b' 
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(a)  Means 

To  make  the  formulas  more  general,  we  consider  here 
that  the  random  fields  x(m,n)  and  T(m,n)  are  not  zero-mean. 
Since  the  fields  are  homogeneous,  the  mean  is  a constant 
and  the  same  for  all  pixels.  All  the  previous  equations  of 
the  image  filters  remain  valid,  provided  we  change  x to 
(x  - x)  and  y to  (y  - x) , where  x is  the  mean  value  of  the 
field  x (m,n) . 

We  will  use  the  following  means: 


E (x)  = x 

E(T)  = T 


E(v)  = 0 


E(y|H0)  = yo  = X 


E (y  1 H1)  = y1  = ax  + T 

A 

E(Xp)  = E(x)  = x 

(b)  Variances 

The  variances  are  given  by 

V (x)  = o2 


V (Xp)  - *p 


(6.5) 
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(c)  Covariances 


We  have  the  following  covariances 


(6 


Cov(v,x)  = Cov(v,x  ) = 0 


Cov(v,y^)  - r , i = 0,1 


Cov(x  ,x)  = b 

ir 


(6 


Cov(Xp,yQ)  = Cov (Xp, x)  = b 


Cov (x  ,y , ) = aCov  (x  , x)  = ab 

p X p 

(d)  Coefficients  of  Correlation 

The  coefficients  of  correlation  are 


Cov(x  , x) 

£ 

v 


v 


Cov(xp,yo) 

Vy, 


a a 

P Y, 
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(6 


Cov(x  ,yL) 


a a 
p yn 


= a 


a a 
p y i 


Now,  returning  to  equation  (6.4)  and  using  the  Gaussian 
pdf ' s , we  have 


L(y,xp) 


[y-E (y  1 xpHQ)  ]2  [y-E(y  IXpHj^)  1 2 


V[y|xpHo] 


V[y|x  H1] 


(6.9) 


where  we  have  dropped  some  constants,  since  they  are  not 

important  in  the  test  against  a threshold. 

/\ 

The  regression  of  the  mean  of  y on  xp  is 


y,- 


E(y|x  Hi]  = y±  + ~ pi(x  - x ) i - 0,1  (6.10) 


Substituting  equations  (6.5)  and  (6.8)  into  (6.10)  gives 


EtylVo1 


= (l-Y)x  + yxr 


E[y|XpH1]  = a (1-y) x+T+aYX^ 


where  we  have  defined: 


(6.11) 


Y = 


a 

°P  Pp 


(6.12) 


Having  the  means  of  the  conditional  pdf's,  we  find 
its  variances: 
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> 


. £ ' -• 


v[y|xpHo] 


Ef [y0-E(y|xpHo) ] } 


ng 

Vis' I Vo1  : 

Using  equations 

A 

Cov  ( X 

fc*' 

A 

v[y|xp 

viylXpH^ 


and 


— ^ 2 

E{ [x+v- (1-y) x-yx  ] > 

ir 

E{ [ (x-x) -y (x  -x) +v] 2 } 

r1 


= a2  + y2a  2 + r - 2yCov(x  ,x) 
P P 

(6.7)  and  (6.12)  we  have 
'*>  = Wp2 

Ho>  ■ °2  - Y0p2  + r 

= El ly1-E(y!xpH1) ]2) 

= E{ [ax+T+v-a (1-y ) x-T-ayx  ] } 

ir 

A 

= E{  [a(x-x)  + (T-T)-ay  (x  -x)+v]' 

b' 


2 2 

a a 


2 2 
a y a. 


v[y|x  hi 


+ a 


+ r 


Substituting  equations  (6.11),  (6.13)  and  (6.14)  into 


(6.9)  gives 


L(Y/Xp) 


_ ~ 2 _ — 2 
[y- (1-y) x-yxp]  [y-a (1-y) x-T-ayxp] 

2 2 2~^  2,  2 2 27“  2 


a -y  a +r 


a(a-ya) +aT  +r 


Let's  define: 


M = (1-y)x 


(6 . If 


2 2 

N = a -YOp 


My,xp) 


~ 2 ^ — 2 
(y-yx  -M)  (y-ayx  -T-aM) 

£ - — £ 


EL 

N+r 


2 2 
a N + r + aT 


(6.1( 


(a2N+r+a  2) (y-yx  -M) 2- (N+r ) (y-ayx  -T-aM) 2 
_£ E 


(N+r) (a2N+r+aT2) 


The  denominator  is  a positive  constant,  because 


N = a2  - y2g  2 = J2(l-p  2)  > 0 

r Jr 


Dropping  the  denominator  and  simplifying  the  numerator 
of  (6.16),  the  likelihood  ratio  becomes 


A /V  ^ ^ A ^ ^ 

L(y,x  ) = [ (y-x) -y (x  -x) ] - 5 [ (y-ax-T) -ay (x  -x) ] 


where 


6 = 


N+r 

2 2 
a N+r+a<r 


(l-pp2)a2  + r 
a2(l-pp2)a  +r+aT2 


(6.18) 


Given  the  prediction  xp  and  the  measurement  y,  the 
decision  rule  consists  in  computing  the  likelihood  ratio 

/v 

Mxp,y)  and  compares  it  against  a threhold  if  greater, 
we  choose  hypothesis  H^,  if  not  we  choose  hypothesis  Hq: 


LUp,y) 


. Hypothesis 


. . . Hypothesis  H 


In  figure  6.3  a block  diagram  of  the  estimator-detector 
filter  is  shown.  It  is  assumed  that  the  background  pre- 
dominates in  the  picture  and  the  objective  is  the  detection 
of  a smaller  texture  that  we  call  target.  The  background 
texture  is  modeled  as  in  chapter  II.  To  accomplish  such 
modeling,  however,  a picture  without  target  is  needed.  In 
an  actual  application  it  is  realistic  to  assume  that 
we  have  a typical  realization  of  this  process  (background) 
before  the  target  enters  in  the  scene.  For 
example,  assume  that  we  are  dealing  with  a surveillance 
problem,  where  some  area  has  been  scanned  many  times  looking 
for  a very  low  contrast  target.  Therefore  a recursive  or 


hybrid  filter  can  be  designed  to  predict  and  estimate  the 
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Estimator-detector  filter 


background.  About  the  target  we  need  only  to  have  a good 
guess  about  the  mean  and  variance  of  its  gray  level.  In 
cases  of  very  low  contrast  the  mean  may  be  assumed  equal 
to  that  of  the  background.  The  variance,  ir  many  cases, 
may  be  assumed  much  smaller  than  that  of  the  background. 

Also,  in  many  cases,  the  observation  noise  is  much  smaller 
than  the  background. 

The  estimator-detector  filter  of  figure  6.3  scans  the 
image  row  by  row  from  top  to  bottom.  Assume  that  the 
previous  pixels  were  correctly  decided  as  belonging  to 
the  background.  The  decision  box  receives  the  measurement 

A 

y and  the  prediction  xp  and  computes  the  likelihood  ratio 

A 

) from  equation  (6.17).  If  this  ratio  is  less  than  a 
chosen  threshold  n , the  pixel  (m,n)  is  considered  background 
and  the  observation  y is  used  to  find  the  background  estimate 
x(m,n).  If  the  ratio  is  greater  than  n,  the  pixel  is  con- 
sidered as  belonging  to  the  target  region.  While  in  the 
target  region  the  background  filter  cannot  use  the  observa- 
tion y,  unless  we  have  "a  priori"  knowledge  about  the  target 
gray  level.  Since  we  are  considering  that  the  target  region 
is  small,  the  background  filter  may  run  in  its  prediction 
mode  without  much  degradation,  provided  the  dynamic  model 
is  accurate  enough.  Observe  that  both  recursive  and  hybrid 
predictors  can  be  used,  the  only  difference  being  the  box 
named  "improved  prediction"  that  applies  only  to  the  hybrid 
predictors.  Since  the  three-dimensional  predictors 
presented  in  chapter  V,  makes  use  of  the  measurements  in 


the  previous  frame,  we  can  expect  a good  performance  of 
those  filters  in  the  target  region. 


* » 

( 


The  gray  level  T(m,n)  of  the  target  texture  can  be 

A 

estimated  using  the  background  prediction  xp  and  the 
observation  y|H^. 


T(m,n)  = k1y1  + k2xp 


where 


y1  = y | H1  = ax  + T + v 


The  estimation  error  is 


*e  = T-T  = T - k1(ax+T+v)  - k2xp  (6.19) 


In  order  to  choose  an  unbiased  estimator: 


E(e)  = T - k^ax+T)  - k2x  = 0 


k-  = — [T-k. (ax+T) ] 

x 


(6.20) 


Using  equations  (6.19)  and  (6.20)  the  weights  k^  and 
k2  can  be  found  using  the  minimum-mean- squared  error  criteria: 


kl  = 


— 2 —2  2—2  2 
axT(ap  -b)  + x cT  +T  Jp  ) 


a2x2 (a2+op2-2b)+2axT(ap2-b)+x2aT2+T2ap2+x2r 


(6.21) 
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Observe  that  if  we  neglect  the  observation  noise 
(r  = 0)  and  the  prediction  error  (x  = x,  b = a 2 = a2) , 
equations  (6.21)  and  (6.20)  become 

k^  = 1 , k2  = -a 


thus : 


T = y1  - axp  (6.22) 

Equation  (6.22)  is  intuitively  appealing,  since  it  is 
exactly  what  we  expect  in  the  absence  of  observation  noise 
and  with  an  ideal  predictor.  The  block  diagram  of  figure 
(6.3)  considers  this  simplified  case. 

D.  SPECIAL  CASES 

In  this  section  we  will  study  some  particular  cases  of 

detection.  The  detector  developed  in  the  previous  section 

will  be  applied  to  simple  cases  of  practical  interest  and 

its  performance  (R.O.C.)  will  be  derived. 

2 

1.  Case  I:  a = 1,  a =0 

In  this  case  the  target  gray  level  is  deterministic 
and  immersed  in  the  background.  The  target  gray  level 
T is  an  unknown  constant  and  we  also  don ' t know  the  target 
shape  or  location  in  the  picture.  The  objective  is  to 
suppress  the  background  in  order  to  enhance  the  target. 

As  before,  the  recursive  or  hubrid  filters  are  used  to 


predict  the  background  x and  the  detector,  using  this 
prediction  and  the  observation  y,  decide  whether  the  pixel 
(m,n)  is  target  or  background.  The  output  is  a binary 
image  where  1^  is  target  and  0_  is  background. 

The  likelihood  ratio,  equation  (6.17),  becomes 

A A A A A 

L(y,x  ) = [ (y-x) -y (x  -x) ] ^ - 5 [ (y-x-T) -y (x  -x) ] 

P P P 


where 


6 = 1 


Thus 


L(y,x  ) = ,2T [ (y-x) -y (x  -x) ] - T^ 

P P 


Dropping  the  constants,  this  ratio  becomes: 


L(y,x  ) = (y-x)  - y (x  -x) 

P P 


(6.23) 


where 


Y a °p 

P 


Equation  (6.23)  is  intuitively  appealing,  because  it  is 


exactly  the  residue  of  the  recursive  filters  when  y = 1 . 


H 


This  would  happen  if  the  filters  were  optimum.  Such  a 
decision  rule  with  y = 1 was  used  in  [5]  as  a generalization 
of  the  one-dimensional  Kalman  filtering  technique.  However, 
we  have  shown  here  that  such  an  extension  is  not  correct, 
because  the  filters  are  not  optimum,  and,  consequently, 
the  residue  (or  innovation  process  [16])  does  not  correspond 
to  the  likelihood  ratio.  However,  it  is  a good  approxima- 
tion to  make  y = 1 for  the  recursive  filters,  in  cases  of 
low  observation  noise  and  highly  correlated  background.  It 
can  be  verified  that  y = 1 for  the  hybrid  filters. 


Now,  we  will  make  the  assumption  that  y and  x have 


a joint  Gaussian  distribution  in  order  to  find  the.  performance 
of  the  detector.  Previously  we  had  assumed  that  the  con- 
ditional pdf's  were  normally  distributed  and  we  validated 
such  assumption  with  experimental  results,  but  this  does 


H 


not  mean  that  the  joint  distribution  is  also  Gaussian. 
However,  it  is  analytically  convenient  here,  in  order  to 
have  simple  results. 


(6.2 


=a  + y a - zyd 
yi  p 


2 2 2 

= a - Y a 

yi  P 


2 2 2 
« a + r - y ap^ 


= (l-pp2)o2  + r 


The  variable  L has  the  same  variance  for  both  hypothesis 
Hq  and  H1,  only  the  means  are  different  as  given  by 


E(L|H0)  = lQ  = (x+v)  - x = 0 


ElLlHj^)  = 4x(x+T+v)  - x = 


In  figure  6.4  the  conditional  pdf's  of  the  likeli 

/A 

hood  ratio  L(y,x  ) is  shown. 

The  probability  of  false  alarm  is 


Pp  = prob (L  j Hq  > n) 


— / exp(-£2/2a  2)dl 


/2tt  n 


— / exp(-u2/2)du 

/2tT  . 

n/aL 
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erf c (n/aL) 


(6.26) 


The  probability  of  detection  is 


pp  = prob(L|H1  > n) 


/ exp[-  (£-T)  2/2oT2]dl 

/Tir  a.  L 

l n 


PD  " 

L 


(6.27) 


where 


erfc (a) 


00 

/ exp(-u^/2)du 

a 


The  choice  of  the  threshold  can  be  done,  for 
example,  by  specifying  the  maximum  admissible  probability 
of  false  alarm;  in  this  case  we  have  a Neyman-Pearson  detector 
[15] . Another  criteria  is  to  minimize  the  overall  probability 
of  error: 


PE  = p(H0)p(ElH0>  + P(H1)P(E|H1) 


P (H  ) Pp  + P(H.  ) (1-PrJ 


If  we  know  P(HQ)  and  P(H1),  an  optimum  threshold 
can  be  found  that  minimizes  Pg.  If  this  is  not  the  case, 
the  best  we  can  do  is  to  assume  that  the  hypotheses  are 
equally  likely,  P(Hq)  = P(H^)  = 1/2.  With  this  choice, 
it  can  be  shown  that  the  optimum  threshold  is  n = T/2  and 

PE  = PF  = 1 PD‘ 

The  implementation  of  the  estimator-detector  for 
this  case  is  similar  to  the  one  shown  in  figure  6.3.  In 
this  case,  since  we  assume  a target  with  constant  gray 
level,  and,  also  assuming  neglectable  observation  noise, 
we  might  use  the  measurement  in  the  recursive  estimation, 

A 

while  the  filter  is  in  the  target  region,  provided  that  T 
is  subtracted  from  yfH^.  The  best  estimate  of  T must  occur 
at  the  edges  of  the  target  region,  because  there,  the 
background  prediction  xp  is  the  best. 

2.  Case  II:  a = 1,  aT2  / 0 

In  this  case  the  target  gray  level  is  a random 

_ 2 

variable  with  known  mean  T and  variance  and,  as  before, 

we  don't  know  its  shape  or  location  in  the  picture. 

The  likelihood  ratio,  equation  (6.17),  becomes 


L(y,x  ) = [ (y-x) -y (x  -x) ] 2 - o [ (y-x) -y (x  -x) -T] 2 


= ( 1—5 ) [ (y-x)-YU  -x)  ] 2 + 26T  [ (y-x) -y  (x  -x)  ] -6T2 

r*  c*' 

(6.28) 


where 


(1-p  2)a2  + r 

2 — 2 1 

(l~Pp  J 


and 


1-6  = 


2 2 2 
(1_pp  )0  +r+0T 


> 0 


Rearranging  equation  (6.28)  in  order  to  have  a 
quadratic  form,  and  observing  that  (1-6)  is  positive: 


L(y>xp) 


(1-5) 


A A -A 

[ (y-x) -Y (x  -x) P + 28 [ (y-x) -y (x  -x) ] 
P P 


+ 8 - 3(3  + T) 


= [(y-x)-y(x  -x)+6]  - 3(3  + T) 

P 


Dropping  the  constants,  this  ratio  becomes 


L(y,x  ) = [ (y-x) -y (x  -x) +8]  (6 

P P 


where 


= T [ ( 1-p  2 ) -%+-%] 


29) 


3 


To  calculate  the  performance  of  the  detector,  we 


make  again  the  assumption  that  y and  x^  have  a joint 


Gaussian  distribution. 

Define  the  random  variable  Z: 


Z = (y-x)  - y(x  -x)  + 3 


Since  y and  xp  have  a bivariate  normal  distribution. 


Z is  also  normally  distributed  and  its  pdf  is 


PZ(z) 


/2tF  Or 


exp  [TJ-Z-zli] 


(6.30) 


2a, 


z.  = y - x + 3 


= °y2  + y2op2  “ 2oi  Y ayap  i = 0,1 


2 2 2 _ 2 2 2 

= a -va  = a -pa 

y p y p 


Since  L = Z , its  pdf  is 


PLU) 


— IP  7(/I)  + P7(-/I) 

2/T  2 2 


Thus 


PLU) 


i {exp[~-^L--zli]  + exp  } 

2a„/2iFT  2az  2aZ 


(6.31) 


/ 


The  mean  and  variance  of  L(y,x  ) can  be  computed 

t' 

using  the  Moment  Generating  Function  of  Z. 


M (t)  = E (ezt ) = exp [ tz  + §^z2t2] 


E (L)  = E ( Z ) = M 


2)  = M 0 ) = z2  + a 2 


V (L)  = E(L2)  - E2(L)  = E(Z4)  - E2(Z2) 


= MZIV(0)  - [MZIZ(0)]2 


Therefore 


E (L)  = z + a, 


V (L)  = 


2a„2(a  2 + 2z2 ) 

u L 


(6.32) 


The  distance  between  the  means  of  L(y,x  ),  under 
the  hypothesis  Hq  and  H^,  is  calculated  using  (6.30)  and 
(6.32)  : 


d = E(L|H1)  - E(L|Hq) 


E(L|H1)  = (T+3) 2+a2+aT2+r+y2ap2-2o1ay^Gp 


— 2 2 2 2 2 
= (T+S)  +o‘+oT  +r-Y  ap 


= (T+S)2+aT2+(l-op2)a2+r 
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E(L|H  ) = 0 + o + r + y o - 2PqY a cp 

* 1 o ' 


= 82  + a2  + r - y2p  2 

fc*' 


o 2 2 

= + (l-p_  )a^  + r 

&r 


thus : 


d = T2  + 2ST  + aT2 


(6. 


From  equation  (6.33)  we  can  see  that  the  distance 
between  the  means  is  non-zero  even  for  a target  with  zero 
mean. 

The  probability  of  detection  is 

00 

PQ  = probUjH^  > n)  = / PL(£|H^)d£ 

n 

Using  equation  (6.31): 


P - i r L 

*^n  2 J 


n z. 


/2irl 


- (/I  - z.  ) 2 

exp  [ 5 ]dl 

2o 


+ 7 / 


-(/I  + z.)2 
exp  [ 5 Ids, 


a /2tTT 


2a 


2 


(6. 


/ 


Using  the  complementary  error  function  (erfc) 
equation  (6.34)  becomes 


PD  = erfc  [^—  (n^2  - z-^) ] + erf c [^— (n1//2  + z1)  ] 

(6.35) 


Z1  Z1 


Similarly,  the  probability  of  false  alarm  is  given 


by 


P„  = ProbUjH  > n) 

F O 


/ PL(*|H0)dU 


(6.36) 


1 


I 

< 


= — (background-to-noise  ratio) 


—2  2 
T +aT 

^ — (target-to-background  ratio) 

a 


Normalizing  the  equations  with  respect  to  the 
2 

background (a  = 1),  equations  (6.35)  and  (6.36)  become 


PD  ~ erfCW_,  TBR 
1+RTT 


/b+RTT] 


(6.39) 


+ erf c [1 


/b+RTT] 


1+RTT 


if- 


/b]  + erfcOM  + /b  ] 


(6.40) 


where 


* - Mr  + (1  ' -V1 


b = a |||(1  + RTT) 


(6.41) 


Using  (6.39)  and  (6.40)  we  have  calculated  the 
R.O.C.  (Receiver  Operating  Characteristic)  of  the  detector. 
In  figure  6.5  the  effect  of  the  target-to-background  ratio 
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(TBR)  is  shown.  As  expected,  the  detector  performance  is 
highly  dependent  of  this  ratio.  For  a probability  of 
false  alarm  of  10%  the  probability  of  detection  varies 
from  67%  to  97%  when  TBR  varies  from  -3  to  3 dB.  In 
figure  (6.6)  the  effect  of  the  background-to-noise  ratio 
is  shown.  Since  this  ratio  is  responsible  for  the  accuracy 

A 

of  the  background  prediction  (x^) , it  is  also  an  important 
factor  in  the  detection.  From  this  figure  we  can  see  that, 
for  a probability  of  false  alarm  of  10%,  the  probability  of 
detection  increased  from  63%  to  93%  when  the  BNR  varied 
from  10  to  30  dB. 

In  figure  (6.7)  the  effect  of  the  randomness  of  the 
target  gray  level,  represented  by  the  ratio  RTT  is  shown. 

We  can  also  observe  a considerable  improvement  in  the 

detection  when  this  ratio  is  increased. 

2 

3.  Case  III:  a = 0 , =0 

In  this  case  the  target  gray  level  is  deterministic 
with  known  value,  though  its  shape  and  location  is  unknown. 
The  target  and  background  regions  are  completely  uncorrelated 
(a  = 0) , thus  it  is  the  case,  for  example,  of  images  in 
the  visible  spectrum  where  there  is  target  or  background, 
but  not  the  addition  of  these  textures . 

The  likelihood  ratio,  equation  (6.17),  for  this 
case,  becomes 

L(y,x  ) = [ (y-x) -Y  (x  -x)  ] 2 - <$[y-T]2  (6.42) 


12  2 
6 = p(o  -yap  +r) 


1 2 2 
±(1-Pp2)a2+1 


Let's  find  the  distance  between  the  means,  under 
the  two  hypotheses  Hq  and  H^: 


A A A A A 

L(y,x  ) = (y-x)  +Y  ix  -x)  -2y (x  -x) (y-x) -6 (y-T) 

P P P 


y H = x + v 
1 o 


E(L|HQ)  = a2+r+y2ap2-2Yb-5E[ (y-x)-(T-x) ] 2 


2 2 2 2 — 2 
= a +r-Y  a -5  (a  +r) -6 (T-x) 

P 


E(L)H  ) = (1-5) (a2+r)  - y2c  2 - 5(T-x)2  (6.43) 

o p 


y | H1  = T + v 


e(l|h1)  = 


(T-x)2  + (1-5) r + '(2°rj2 

P 


(6.44) 


Subtracting  (6.43)  from  (6.44)  we  find  the  distance 


d = E(L|HX)  - E (L | Hq) 


d = (1+6) (T-x)2  + ( 2p  2+6-l) c2 

P 


(6.45) 


Observe  that  the  distance  is  non-zero,  even  for  tar- 
gets with  the  same  gray  level  as  the  mean  of  the  background 
(T  = x)  . 


This  case  might  cover  the  model  of  some  images  with 


two  quite  distinct  textures.  A difficult  problem,  for 
example,  is  the  restoration  of  pictures  with  two  different 
regions,  one  is  the  sky  and  the  other  the  sea.  This  is  a 
typical  case  where  we  cannot  model  the  image  as  a homogeneous 
random  field.  One  approach  to  this  problem  is  to  model  each 
texture,  separately,  as  two  homogeneous  random  fields,  by 
the  method  of  Chapter  II,  and  design  one  image  filter  for 
each  texture.  To  apply  these  filters,  first  we  have  to 
decide  to  which  texture  the  pixel,  to  be  processed,  belongs. 
Such  a decision  can  be  made  by  the  detector  of  the  present 
case.  In  figure  (6.8)  the  block  diagram  for  this  situation 
is  shown.  Following  the  same  notation  as  before,  we  have 
a filter  for  the  background  x (one  texture)  and  another 
filter  for  the  target  T (the  other  texture) . The  threshold 
decision  selects  the  proper  filter,  while  the  other  filter 
runs  in  its  prediction  mode. 

The  likelihood  ratio  for  this  case,  equation  (6.17), 
is 


My,x  ) 


[ (y-x) -y (x  -x) ] 2 - 6 (y-T) 2 


(6.46) 


Let's  compute  the  distance  between  the  means: 


E(L|H  ) = (1-6) (c2+r)-y2a  2-5 (T-x) 2 

ir 

E(L|H1)  = (1-6) (cT2+r)  + y2ap2  + (T-x) 2 

d = E(L|H1)  - E(L|HQ) 

d = ( 1+6 ) (T  - x)  2 + (1-6)  (a  2-a)  +2p  2a2  (6.47) 

Equation  (6.47)  shows  clearly  the  differences 
between  the  two  textures  (mean  and  variance) . 

As  in  Case  IIT,  the  performance  of  this  detector 


can  be  found  using  the  Gaussian  hypothesis,  however  it  is 
not  possible  to  have  an  elegant  analytical  solution. 


Figure  6.8 

Estimator-detector  filter  for  two  textures 


VII.  TARGET  TRACKING 


A.  INTRODUCTION 

In  the  previous  chapters  we  have  developed  statistical 
image  filters  that  exploit  spatial  and/or  temporal  redun- 
dancies existing  in  the  pictures.  One  application  of  those 
filters  is  the  restoration  of  images  contaminated  with 
observation  noise.  In  this  application,  the  objective  is 
to  remove  the  observation  noise,  in  order  to  recover  the 
original  image.  Another  application  is  to  extract  some 
desired  information  contained  in  the  pictures,  which  we  have 
called  target.  In  this  case,  we  want  to  suppress  the  rest 
of  the  information,  which  we  have  called  background,  in 
order  to  enhance  the  target. 

The  target  detection  and  tracking  problem,  by  means 
of  pictures  sequenced  in  time,  may  be  divided  into  three 
phases : 

(a)  Extraction  of  the  targets  from  the  background, 
creating  a binary  picture  (1-target,  O-background) . 

(b)  Recognition  of  the  target (s)  that  we  are  interested 
in  tracking. 

(c)  Track  the  target (s)  from  frame  to  frame. 

This  research  addressed  the  first  and  third  phases. 

The  recognition  problem,  for  our  present  purpose,  will  be 
considered  as  solved  by  a human  operator,  or,  in  the  case 
of  an  automatic  system,  by  some  existing  algorithms  [18J , 
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[19]  of  pattern  recognition.  The  extraction  problem  was 
covered  by  the  previous  chapters. 

In  this  chapter  we  address  the  tracking  of  targets 
from  frame  to  frame.  The  basic  idea  is  to  model  the  target 
movement  within  the  pictures,  and  to  construct  a Kalman 
filter  to  track  the  target  centroid,  using  the  image  filters 
as  the  measurement  device.  Of  course,  others  relevant 
points  of  the  target  might  be  tracked,  in  a similar  way, 
in  order  to  find,  for  example,  the  target  attitude.  The 
image  filters  referred  to  here  are  the  three-dimensional 
recursive  and  hybrid  predictors,  presented  in  Chapters  IV 
and  V,  respectively,  and  the  detector  developed  in  Chapter 
VI. 


B.  TARGET  TRACKING  FILTER 

In  this  section  we  will  construct  a Kalman  filter  to 
track  the  target  centroid.  Assume  that  we  receive  pictures 
at  some  specified  rate,  say  30  pictures/sec.  Since  the 
pictures  have  only  two  dimensions,  the  target-image  dynamics 
is  related  only  to  the  components  of  the  actual  movement 
that  are  parallel  to  the  picture-plane  (see  figure  7.1). 

The  effect  of  the  movement  perpendicular  to  the  picture  is 
.to  increase  or  decrease  the  size  of  the  target-image.  If 
the  target  has  the  same  aspect  at  all  angles  (for  example: 
a sphere) , then  it  is  possible  to  compute  such  component  , 
using  successive  frames. 
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The  target-image  becomes  a point-target  at  large  dis- 
tances, but,  in  general,  it  is  a mass-target  with  some  shape 
which  changes  when  it  maneuvers . Such  change  can  have  very 


interesting  applications  in  the  target  tracking  problem, 
since  it  gives  useful  information  of  target  maneuvering. 

In  what  follows,  we  model  the  movement  of  the  target 
centroid  in  the  plane  of  the  picture  and  construct  a Kalman 
filter  to  estimate  its  x and  y coordinates,  and  also  to 
predict  its  position  in  the  next  frame. 


> 


1.  Dynamic  Model 

Since  we  don't  know  the  intention  of  the  target 
driver,  we  will  consider  two  random  accelerations  in  the  x 
and  y directions,  therefore  the  position  of  the  target- 
image  is  given  by 


x(k+l)  = x(k)  + x(k)T  + ^w^(k)T2 


y(k+l)  = y(k)  + y(k)T  + ^w2(k)T2 


(7.1) 


Equation  (7.1)  assumes  that  the  time  separation  T 
’between  frames  is  small  enough  such  that  we  can  consider 
that  the  acceleration  and  velocities  are  constants  between 


frames . 


The  velocities  are  given  by 


x(k+l)  = x(k)  + w^(k)1 


y (k+1)  = y (k)  + w2  (k)  T 


(7.2) 


where  we  have  used  in  this  derivation  the  simplified 
notation  x(k)  for  x(kT). 

To  use  state  variable  notation,  define  the  state 
vector  x and  the  random  forcing  input  w: 
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The  centroid  of  the  target-image  is  measured  by 
the  image  filters  previously  developed.  The  pixels  are 
examined  one  by  one  and  classified  as  target  or  background. 
After  this  classification  the  target  centroid  is  computed 
by 
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xm  * I l nkT(mk'nk>  ' s * l T(VV 
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(7.5) 


where  (m^,n^)  are  the  coordinates  of  the  pixels  classified 
as  target.  If  we  create  a binary  picture  with  T(mk,n^)  = 1, 
equations  (7.5)  become 
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(7.6) 
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Since  the  image  filters  are  not  perfect,  some 
target-pixels  are  missed  and  some  background-pixels  are 
taken  as  targets.  Therefore  the  measurement  of  the  centroid 
has  an  error: 


is: 


(7.7) 


y = y + v_ 

* 2 

where  v1  and  v2  are  the  errors  in  the  x and  y coordinates, 
respectively. 

Let's  analyze  the  errors  v1  and  v2 . First  deter-- 
mine  if  they  are  uncorrelated.  In  figure  7.2  assume  that 
the  image  filter  missed  the  target-pixels  in  region  A. 

The  result  is  an  error  in  the  observation  of  the  centroid. 

Both  x and  y coordinates  will  be  in  error.  In  other  situa- 
tions only  one  coordinate  is  wrong.  In  general,  we  can  not 
draw  any  inference  about  one  error  by  the  knowledge  of  the 
other,  therefore  we  will  assume  that  v^  and  v2  are  uncorrelated 


Figure  7.2 


Measurement  error 


We  must  determine  if  the  errors  are  correlated  with 
the  coordinates  (x,y) . Since  the  background  and  the  target 
gray  levels  are  considered  homogeneous  random  fields, 
there  is  no  reason  to  think  that  the  image  filter  will  be 
more  likely  to  make  errors  for  some  positions  of  the 
target  within  the  frame,  or  at  some  special  frames.  On 
the  other  hand,  the  nature  of  the  filter  itself  could  be 
a source  of  correlation.  Both  image  filters  (recursive 
and  hybrid)  use  recursion  in  space,  therefore  the  background 
prediction  is  poor  near  the  top  left  corner  and  also  near 
the  first  row  and  column  of  all  frames,  because  of  the 
boundary  transients  of  the  filters.  In  the  first  frame  the 
hybrid  predictor  does  not  make  use  of  the  time  correlation 
in  this  frame.  In  the  recursive  predictor  the  first  few 
frames  have  transients  in  recursion  in  time.  The  conclusion 
is  that  there  are  some  positions  (x,y)  in  the  time-frame 
that  have  poor  probability  of  detection.  However,  we  have 
seen  that  the  recursive  filters  converge  very  fast,  there- 
fore, we  will  neglect  such  correlation,  since  it  occurs  only 
in  a very  small  region. 

One  must  test  to  see  if  the  errors  are  zero-mean. 

The  detection  error  could  happen  at  any  pixel,  but  the  bias 
is  most  likely  in  the  direction  that  the  picture  is  scanned 
(top  left  corner  to  right  bottom  corner) , because  when  the 
filter  is  in  the  target  region  the  background  predictors, 
in  general,  can  not  use  the  observation  y(m,n),  which 
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results  in  a degradation  in  the  detection.  However,  we 
expect  that  such  degradation  can  be  neglected  for  practical 
applications  where  the  target  has  small  size  and  the  back- 
ground correlation  is  high.  There  are  other  schemes  for  pro- 
cessing that  can  be  done  to  overcome  this  problem.  One  is 
for  example,  to  process  the  pictures  top-to-bottom  and 
bottom-to-top,  alternatively,  from  frame  to  frame. 

After  these  considerations  we  conclude  that  it  is 
reasonable  to  consider  the  errors  v^  and  v2  zero-mean  white 
noise  and  uncorrelated  with  x and  y.  Such  conclusions 
simplify  considerably  the  following  derivation  and  the 
result  will  be  a very  simple  Kalman  filter. 

Define  the  measurement  and  error  vectors  as  below 
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In  state  space  notation,  equation  (7.7)  becomes 


£(k)  = 


x (k)  + v (k) 


(7.9) 


where  v(k)  is  a zero-mean  stationary  random  process,  uncor- 
related with  x(k)  and  with  the  covariance  matrix  below 
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The  dynamic  model  is 
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The  measurement  model  is  given  by 


£(k)  = H x(k)  + v (k) 


(7.15) 


where 


(7.16) 


(7.17) 


The  random  forcing  input  w(k)  and  the  measurement 
noise  v(k)  are  considered  white  noise  uncorrelated  between 
themselves  and  also  with  the  states.  Its  covariance  matrices 
are  given  by 
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(7.18) 


The  Kalman  filter  equations  are  [17]  : 


x (k)  = x' (k)  +G(k)[£(k)  - Hx' (k) ] (7.19) 
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where  x(k)  is  the  estimated  state  vector  and  x' (k)  is  the 
predicted  state  vector  given  by 


x' (k)  = $ x (k-1) 


(7  r 20 ) 


G(k)  is  the  filter  gain,  computed  as  below  [17] : 


G (k)  = P’  (k)  HT[HP'  (kJH^R]"1 


(7.21) 


P (k)  = [I-G  (k)  H]  P ' (k) 


(7.22) 


p' (k+i)  = $p(k) $T  + rorT 


(7.23) 


where  P(k)  is  the  covariance  matrix  of  the  state  error 
and  P' (k)  is  the  predicted  covariance  matrix  of  the  state 
error. 

The  estimated  state  error  vector  is 


e(k)  = x(k)  - x(k) 


The  covariance  matrix  of  the  state  error  is 


P(k)  = E [e (k) e (k) ] 


The  predicted  state  error  vector  is 


(k)  = x'  (k)  - x (k) 


(7.24) 


(7.25) 


(7.26) 
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The  predicted  covariance  matrix  of  the  state  error  is 


P'(k)  = E [e ' (k)  e ' T (k)  ] 


(7.27) 


The  gain  G(k)  is  computed  by  initialization  of 
equation  (7.21)  with  P'(l)  to  have  G(l),  then  find  P(l) 
by  equation  (7.22)  and  P'(2)  by  equation  (7.23).  Following 
this  procedure  the  matrices  G(k),  P(k)  and  P' (k)  are 
calculated. 

Since  the  x and  y coordinates  of  the  target  centroid, 
as  well  as  its  measurement  errors  v^  and  are  uncorre- 
lated, we  can  initialize  equation  (7.21)  with  the  predicted 
covariance  of  the  state  error  below 
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matrices  G(k) 

and 
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Pn(k) 

tl-gi(k)  lp*_1(k) 

P12  ^ 

= 

[l-g-L  (k)  ] p^2  (k) 

p22  ^ 

= 

p22 (k)  - g2(k)P|2(k) 

(7.32) 

p33(k) 

= 

[l-g3(k)]P'3(k) 

P34(k) 
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P44(k) 
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Initializing  the  filter  equation  with  the  initial 
prediction  x' (1) , we  have  the  estimated  state  vector 


x(k)  = x' (k)  + gL(k) [xm(k)-x' (k) ] 


x (k)  = x' (k)  + g, (k) [x  (k)-x' (k) ] 

^ m 


(7.34) 


y (k)  = y'(k)  + g3 (k) [ym(k) -y ' (k) ] 


y (k)  = y ' (k)  + g4 (k) [ym(k) -y ' (k) ] 


The  predicted  state  vector  is 


x' (k+1)  = x(k)  + T x(k) 


x' (k+1)  = x(k) 


(7.35) 


y' (k+1)  = y(k)  + T y(k) 


y ' (k+1)  = y (k) 


In  figure  7.3  we  show  a flowchart  for  implementation 


of  the  filter. 
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Figure  7.3 

Kalman  filter  flowchart 


C.  FILTER  APPLICATION 


The  tracking  filter  developed  in  this  chapter  is  used, 
primarily,  to  provide  the  sensor  servo  with  the  predicted 
centroid  coordinates  for  the  next  frame.  Besides  this 
application,  the  information  of  the  predicted  position  of 
the  target  can  be  used: 

(a)  To  reduce  the  image  processing, 

(b)  To  improve  the  performance  of  the  target  detector. 

The  reduction  in  the  image  processing  can  be  achieved 

by  delimiting  to  a small  area  to  be  processed,  centered 
at  the  predicted  position  of  the  target  centroid.  Instead 
of  processing,  for  example,  a 100  x 100  picture  we  might 
work  only  in  a 20  x 20  region,  depending  on  the  target  size. 
With  such  reduction  we  could  afford  to  make  more  processing 
in  this  smaller  area,  in  order  to  improve  the  performance 
of  the  target  detector. 

Some  examples  of  this  extra  processing  are: 

(a)  Scan  the  picture  in  two  directions,  one  starting 

at  the  top  left  corner  and  the  other  at  the  bottom 
right  corner.  The  results  are  two  binary  pictures 
(1-target,  0-background)  with  greater  probability 
of  error  for  the  pixels  located  in  the  transient 
region  of  the  recursion  in  space  (see  figure  7.4). 

One  simple  way  to  combine  the  two  estimates  is  to 
pick  the  half  of  each  picture  opposite  to  the  starting 
point,  because,  in  this  way  we  remove  most  of  the 
unfavorable  region  of  detection. 
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If  more  processing  is  allowable  one  might  scan  the 
picture  four  times  starting  at  each  corner  and  pick 
the  quarter  of  each  estimated  picture  opposite  to  the 
starting  point.  In  this  case  the  detection  is 
completely  free  of  the  transient  region. 

(b)  In  the  detection  of  the  target  edges  the  direction 
of  scanning  the  picture  is  very  important.  Since  a 
good  prediction  of  the  background  is  needed  and,  in 
general,  such  prediction  is  degraded  while  in  the 
target  region,  the  result  is  that  the  detection  of  the 
edge  at  the  left  is  better  than  the  one  at  right 
(see  figure  7.5).  One  way  to  overcome  this  problem 
is  to  use  the  estimator-detector  only  in  favorable 
conditions.  Assume  that  the  prediction  of  the  target 
centroid  is  good.  Process  first  the  left  hand  part 
of  the  picture  until  the  line  AB  (see  figure  7.5). 

Next,  process  the  right  hand  part,  starting  at  the 
top  right  corner.  In  this  way,  the  filter  .vill  work 
in  more  favorable  conditions,  since  it  will  face  both 
left  and  right  edges  with  good  background  predictions. 
The  implementation  of  such  a scheme  may  be  simplified 
if  we  invert  the  columns  of  the  right  hand  part  before 
and  after  the  processing.  In  this  way,  the  filter 
estimates  the  left  hand  part  and  it  is  reinitialized 
at  the  boundary  (line  AB)  to  proceed  to  the  estimation 
of  the  inverted  right  hand  part.  Therefore,  the  only 
additional  process  would  be  the  required  inversion. 


VIII.  COMPUTER  SIMULATION 

A.  DESCRIPTION  OF  THE  SYSTEM 

A computer  simulation  of  the  target  tracking  and  detection 
problem  was  implemented.  The  block  diagram  of  the  system  is 
presented  in  figure  8.1.  It  is  composed  of  three  parts: 
time-frame  generation,  image  processing  and  target  tracking. 
The  generator  of  the  time-frame  creates  images  at  some  speci- 
fied rate  and  characteristics.  The  images  contain  the  target 
and  are  contaminated  with  additive  white  Gaussian  noise.  The 
target  has  a random  movement  from  frame  to  frame.  The  images 
are  processed  in  two  modes:  search  and  window.  The  search 
mode  is  used  to  acquire  the  target.  In  this  mode  the  whole 
picture  is  processed  in  order  to  extract  the  target  from  the 
background.  It  is  assumed  that  a recognition  phase  exists 
to  select  the  desired  target.  In  this  simulation  we  have 
used  only  one  target  to  avoid  such  a phase.  After  the 
acquisition  phase,  the  images  are  processed  only  in  a window 
centered  at  the  predicted  position  of  the  target  centroid. 

The  output  of  the  image  filter  is  the  measurement  of  the 
coordinates  Cam/8m)  of  the  target  centroid.  The  target 
tracking  filter  receives  these  measurements  and  updates 

/v  /v 

the  estimation  of  the  target  centroid  (a,S)  and  also 
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Block  diagram  of  the  simulation 
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computes  the  prediction  for  the  next  frame  (a Bp)  which 
is  feedback  to  the  image  filter. 

More  details  about  each  block  of  figure  8.1  are  given 
in  the  rest  of  this  section. 

1.  Time-frame  Generation 

The  background  has  the  autocorrelation  function: 


R(i,  j ,k)  = a2  pj1!  ph  I ^ I pjkl  i,j,k  = 0,±1,±2, 


The  dynamic  model  that  generates  a random  field 
with  such  an  autocorrelation  function  is 


x(m,n,t)  = A x + w(m,n,t) 


where : 


A = [Pv  PhPt  -PvPh  -PvPt  -PhPt  PvPhPtl 


x = [x(m-l,n,t)  x(m,n-l,t)  x(m,n,t-l)  x(m-l,n-l,t) 


x(m-l,n,t-l)  x(m,n-l,t-l)  x(m-l,n-l,t-l) ] 


E [w2 (m,n, t) ] = Q = a2 (1  - py2)  (1  - ph2)  (1  - pt2) 


The  block  "Background  generator"  implements  this 
dynamic  model  by  driving  the  equation  with  zero-mean  white 
Gaussian  noise  and  variance  Q. 
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The  block  diagram  of  the  combined  estimator-detector 
filter  is  presented  in  figure  6.3,  particularized  to  case 
I (a  = 1,  aT2  = 0) . 

The  constant  y was  defined  (see  equations  6.7 
and  6.12)  by 


Y 


COV (x,xp) 


p 


p 


where : 


°p2  = E[V]  = E[Vx+ep)] 


= cov(x,x  ) + E[x  e ] 


For  the  recursive  filter  the  prediction  error  e 


is  not  orthogonal  to  xp,  because  the  filter  is  not  optimum. 


However,  the  hybrid  filter  is  such  that  the  prediction  error 


e is  orthogonal  to  x , therefore: 

sr  c 


y 


for  the  hybrid  filter 


for  the  recursive  filter 
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The  threshold  n is  chosen  such  that  the  probability 
of  a false  alarm  is  below  some  desired  value.  From  equation 
(6.26) 


Pp  = erfc  (J-) 
F aL 


For  n = 5aT  Pp  = 1.5x10 


-6 


In  the  simulation  the  threshold  can  be  specified  by 
the  operator  or  automatically  set  at  the  previous  value 

n = 5ol  = 5 j/(l-  Pp2)  a2  + r (see  equation  6.25) 

Observe  that  this  criteria  is  completely  independent 
of  "a  priori"  knowledge  of  the  target  gray  level. 

We  have  defined  the  following  ratios: 

TBR^  = Taryet-to-background  ratio  at  input 

TBRq  = Target-to-background  ratio  at  output 

BNR  = Background-to-noise  ratio 

TBR 

o 


P.G. 


TUP 


Processing  gain 


P.G.  = 


a2  + r 


BNR  + 1 

(1  - Pp2)  BNR  + 1 


(8.3) 


From  (8.3)  we  can  see  that  one  upper  bound  for 


P.G.  is 


P.G.  < BNR  + 1 


Therefore,  the  observation  noise  has  to  be  well 
below  the  background  in  order  to  have  a substantial  target 
enhancement.  However,  in  many  applications,  the  observation 
noise  is  really  neglectable  when  compared  with  the  back- 
ground and,  therefore,  substantial  target  enhancement  (or 
background  suppression)  can  be  achieved  with  such  a method. 

After  the  detection  of  each  target-pixel,  the 
centroid  is  computed  by 


_ 1 r 

N £ ni 


= £ l rn. 

N L x 


where  (m^n^)  are  the  coordinates  of  the  i target-pixel, 
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. Target  Tracking  Filter 


The  target  tracking  filter  is  the  implementation  of 

the  Kalman  filter  constructed  in  Chapter  VII.  It  receives 

the  target  centroid  measurement  (a  ,8  ) from  the  image  filter 

mm3 

and  computes  the  estimated  position  (a,  8),  as  well  as  the 

predicted  position  in  the  next  frame  (a  ,8  ).  The  filter 

P P 

gains  are  computed  on-line  according  to  the  flowchart  of 
figure  7.3. 

B . SELECTED  RESULTS 

In  this  section  we  present  some  relevant  results  of 
the  simulation.  The  pictures  have  the  size  34  x 70,  compati- 
ble with  the  Tektronics  4012  Display.  The  images  are 
pictorially  displayed  with  8 levels  of  quantization  obtained 
by  superimposing  some  characters.  The  target  has  a 
"diamond"  shape,  but  it  can  be  easily  changed  to  other 
shapes  enclosed  by  a 7x7  matrix.  In  the  tracking  mode 
the  image  processing  is  only  applied  to  a 20  x 20  window 
centered  at  the  predicted  centroid  of  the  target. 

1.  Background  Prediction 

Since  the  decision  rule  is  highly  dependent  on  the 
background  prediction,  we  have  measured  the  prediction 
error  for  several  situations  and  compared  with  the  theoreti- 
cal results.  The  results  are  presented  in  tables  8.1 
through  8.3.  Observe  that  the  theoretical  and  experimental 
results  have  the  same  order  of  magnitude.  We  can  also 
observe  that  the  experimental  results  are  always  greater 
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TABLE  8.1 


NOISE 

VARIANCE 

3D- RECURSIVE 

PREDICTOR 

3D-HYBRID 

PREDICTOR 

THEORY 

EXPERIMENT 

THEORY 

EXPERIMENT 

o 

• 

H 

.293 

.306 

.282 

.360 

00 

• 

o 

.267 

.280 

.258 

.332 

0.6 

.236 

.255 

.229 

.308 

0.4 

.197 

.216 

.195 

.266 

0.2 

.144 

.166 

.150 

. 206 

PICTURE  VARIANCE  = 1.0 


pt  = 0.95 


TABLE  8.2 


NOISE 

VARIANCE 

3D- RECURSIVE 

PREDICTOR 

— | 

3D- HYBRID 

PREDICTOR 

THEORY 

EXPERIMENT 

THEORY 

EXPERIMENT 

1.0 

.351 

.364 

.303 

.353 

00 

• 

o 

.324 

.343 

.283 

.329 

0.6 

.293 

.316 

.260 

.295 

0.4 

.251 

. 266 

.231 

.269 

0.2 

.191 

.207 

.192 

.217 

PICTURE  VARIANCE  = 1.0 

pv  - ph  = °'8 


0.8 


i 


TABLE  8.3 


NOISE 

VARIANCE 

3D- RECURSIVE 

PREDICTOR 

3D-HYBRID 

PREDICTOR 

THEORY 

EXPERIMENT 

THEORY 

EXPERIMENT 

1.0 

.376 

. 388 

. 338 

.378 

0.8 

.349 

.353 

.317 

.352 

0.6 

.317 

. 337 

.292 

.330 

0.4 

.274 

. 294 

.260 

.295 

0.  2 

.213 

.231 

.216 

.244 

PICTURE  VARIANCE  = 1.0 


pv  " ph  = °*8 


Pt  = 0.7 


than  the  theoretical  ones.  Such  discrepancy  is  due  to  the 
fact  that  the  theoretical  values  are  computed  for  the 
steady-state  condition,  but  the  experimental  values  are 
the  result  of  averaging  over  a small  ensemble  (20  frames) 
and  small  images  (34  x 70) , which  means  a reasonable  influence 
of  the  transient  of  the  image  filters. 

Therefore,  we  consider  that  these  results  are  a 
reasonable  validation  for  the  theoretical  methods  of  filter 
design  and  performance  evaluation  presented  before. 

2 . Transition  Region 

To  observe  the  transient  of  the  image  filters,  we 
have  applied  the  decision  rule  to  a time  sequence  of  images 


i 
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tv 


(time-frame)  without  target.  The  threshold  was  set  at 
n = 0.35.  The  characteristics  of  the  images  are: 


a =1,  Pv  = Ph  = 0.93,  p = 0.95,  BNR  = 30  dB 


Li 


The  results  for  the  recursive  filter  are  shown  in  figures 
8.2  through  8.5.  As  mentioned  before,  the  3D-recursive 
filter  gives  degraded  predictions  at  the  first  few  frames 
and  near  the  first  row  and  column  for  all  frames,  due  to 
the  time  and  space  recursion,  respectively.  Observe  that 
in  frame  2 (figure  8.2)  the  false  alarms  occurred  mostly 
near  the  first  column  and  row,  but  they  can  also  be  seen 
in  the  rest  of  the  picture.  In  frame  10  (figure  8.4)  the 
situation  is  better,  because  it  the  steady-state  was  reached 
in  time,  but  the  false  alarms  persist  due  to  the  transient 
in  space.  In  chapter  VII  we  have  suggested  some  methods 
to  reduce  such  undesirable  effects  of  the  transition  region. 

In  this  simulation  we  have  used  a simple  method  that  is 
to  increase  the  threshold  for  the  pixels  at  the  first  5 
rows  and  columns.  The  results  of  such  a method  are  shown 
in  figures  8.3  and  8.5.  Observe  the  substantial  improvement 
in  both  frames  2 and  10.  The  false  alarms  of  frame  2 
(figure  8.3)  are  basically  due  to  the  transient  in  time. 

The  percentage  of  false  alarm  after  frame  5,  as  in  frame 
10  (figure  8.5),  is  quite  close  to  the  theoretical  value. 

Of  course,  the  probability  of  detection  is  reduced  in  the 
transient  region,  due  to  the  higher  threhold,  however  this 
region  is  very  small  for  the  usual  size  of  pictures  (100  x 100) 


False  alarms  of  the  recursive 
filter  at  frame  2,  using 
two  thresholds. 


Figure  8.4  False  alarms  of  the  recursive 
filter  at  frame  10,  using 
one  threshold. 


Figure  8.5  False  alarms  of  the  recursive 
filter  at  frame  10,  using 
two  thresholds. 
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The  results  for  the  hybrid  filter  are  shown  in 
figures  8.6  through  8.9.  Since  this  filter  does  not  use 
recursion  in  time,  its  transient  is  only  in  space.  Simi- 
lar improvement  was  achieved  by  using  two  threholds,  as 
can  be  seen  in  figures  8.7  and  8.9,  compared  against  figures 
8.6  and  8.8.  Observe  that  frames  2 and  10  (figures  8.7 
and  8.9)  have  almost  the  same  number  of  false  alarms,  as 
should  be  expected,  since  there  is  no  transient  in  time. 

Comparing  these  results  we  see  that  the  recursive 
filter  presented  better  performance  than  the  hybrid  filter. 
This  result  is  according  to  the  theoretical  and  experimental 
values  of  tables  8.1  through  8.3,  where  the  recursive  filter 
has  better  performance  for  high  background-to-noise  ratios, 
as  is  the  case  (BNR  = 30  dB) . 

3 . Background  Suppression 

In  what  follows  we  present  some  examples  of  target 
enhancement  (or  background  suppression) . Several  parameters 
can  be  varied  to  evaluate  the  estimator-detectors  performance. 
In  these  examples  we  have  used  the  same  kind  of  background, 
but  different  values  of  target-to-background  ratios. 

The  characteristics  of  the  background  are 

a2  =1,  Pv  = ph  = 0.93,  Pt  = 0.95,  BNR  = 30  dB 

The  correlations  in  space  (Pv,Pjj)  have  the  same  order  of 
magnitude  as  those  observed  in  chapter  II  for  several  real- 
life  pictures.  The  background-to-noise  ratio  (BNR)  was 
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Figure  8.8 


False  alarms  of  the  hybrid 
filter  at  frame  10,  using 
one  threshold. 


Figure  8.9  False  alarms  of  the  hybrid 
filter  at  frame  10,  using 
two  thresholds . 
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chosen  so  high,  because  we  are  interested  here  in  those 
situations  where  the  observation  noise  is  neglectable  when 
compared  with  the  background,  which  is  the  undesirable 
texture  to  be  suppressed,  in  order  o enhance  the  target. 

The  performance  of  the  filters  will  be  given  by 
the  ratios  defined  in  equations  (8.1)  through  (8.3). 

Example  1 

Filter:  Recursive 

TBRi  = 0 dB 

Threshold  = 0.5 

Figures:  8.10  through  8.15 

In  table  8.4  the  performance  of  the  filter  at 
several  frames  is  shown.  Observe  that  at  frame  1 the  3D- 
recursive  filter  reduces  to  a 2D-recursive  filter  which  does 
not  use  the  correlation  in  time.  Comparing  figure  8.11 
with  8 . 13  and  8.15  we  can  see  the  great  improvement  resulting 
from  the  exploitation  of  the  correlation  in  time. 

Example  2 

Filter:  Hybrid 

TBRi  = 0 dB 

Threshold  = 0.5 


Figures:  8.16  and  8.17 
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Figure  8.11  Frame  1 after  the  recursive  filter 
(TBRq  = 12  dB) 


LU. 

LLLLvl33miQ3g333a33333L333l[-Ll.L33L3333ll3l.t-31.\sLL3333L3003g00T0a0aBO 
3H.LL33iaaaaa3333333333LLL3LLLLH.3330Un0333333U.LLL33333330333330 30000 
v ■ NLL33L3333LLLL33LLLM.LU.LLLLLL333B3LL3LLULWLLLU33333030000C3CXXX30 
VN  N SVLLLLLLL3LLL333LU.NLLLNNLXSV  -L333L\'LL'vLss  w(.H.333333333330300aO 
LL\NS\LLLLNLLL3LLL333L3LL33LLLLL^LL3303LNLLLL'VLLU3atJlJOao:X>OaCBaOaXI 

ULn^M.LL3LLL333LLL33333LL3L^LLL^  'v  ^LL333LLLLLLv'>'>LL330Q300Q3Q(XKBQQQ30 
LL'vNv-^LL3LLL333333333333333LLL3t.LLLL33D3I.LLLLL^VNvL1.30333333OOOa3UUOgO 

u.\ 

L3L 

3LL 

\ 

\ 

\ 

M. 

M. 

L 

3 

L 

L 

3 

L 

3 

1 

3 

L JJZZJ 

3 man 

33303 

aajjJunmujJJLLLJL 

annnnnnnnnnini  t * n 

Unnil3n3D3X333LL333 

ULLJDm 

LLL3Q0I 

L3L300I 

13 

13 

13 

3L3JLLNSNM.J 

3L33L3LV-LL30 

3L33ULLNLL3D 

• :,ljB)lp^jB'«lB:B]B  11133,1 

a 3 r|t|ajBjap|B>t|B|B  3.»f«  a 3 B 
i:r  i riB|i|i  = :Tz;33333 

33\ 

\\ 

u 

L 

3 

L 

3 

a 

33303 

5 

• 

33L 

L 

LL 

3 

TTT 

TtP 

!! 

t 

BlBiBlB  Z[B  - r*  1«1* 

331 

2P 

» 

II  1 

k i » 

T 

• ^ : UlitH 

33L 

33L 

p 

| 

g 

3 

rt 

I- 

P 

ip 

I-L-  • » • ■ • i. 

¥■ 

EH 

i i i 

r 

i 

: 

n 

• MB  B 

1111 

■ 

t 

; 

1 

g 

• 

iT: 

*T: 

t T-  f:  1- f 

i:T:j 

[:]»  i)i|i  b[i  i » i 

i * 

V 

a 

Bill 

It 

n 

El 

n 

El 

El 

Cl 

| 

i 

s 

| 

t* 

T:  • 

Sir 

s 

[:T«  »i»r* 

jiti  feh  4<*  ■ ■ 

£| 

rtr 

ip » 

B 

• 

i 

i 

M 

Bill 

till 

ftrfc 

2 

E 

2 

! 

L 

: 

: i 

TTrrTTTnTTTT 

FI*  tf  t 1 * * 

E " 

i > i 

1 

B 

lU 

B B I B 

S 

? 

| 

? 

■ 

f 

E 

r 

IS 

i»t»  »ti  ip  • • * 

fj.crctt  i * i 

s 

*7 

■ 

■ 

• 

5 

Mil 

• I B 1 

P 

P 

: 

? 

? r 

" : 

-I- 1- l:  1*  1 * ! 

■mrnrfr? 

sE " 

»:  r 

B 

I 

ttl 

iiil 

ill; 

" 

r 

" 

* 

? 

ill. 

trmrrv 

rtr  r 

•t»  II 

T 

i 

frtj 

Mil 

:n|:1 

u 

3 

T 

TtT 

: f:  f.  J»  • 

‘ L 

• 

■ 

r 

jakj 

1 1 I 1 

Si 

" 

_ 

"TT 

w8t 

ST 

•R  T 

•1?  » 

t 

B 

ttl 

1 1 1 1 

JTl 

5 

? 

: : 

: :F:TST  » • ST: 

!r  ' 

n?  i 

• 

1 

ttl 

•p»  1 

itrT? 

T 

T 

" " ; 

1 " ; 

1-1-  l-l-l- 

LjPT 

•i»  • 

■ 

a 

rp 

1 | | p 

- 

3 

n 

“ 

TT 

prrfT 

FS 

if:  i i|ij;ji  « > if1 

T? 

rrr 

i 

? 

« 1 1 1 

TT 

£ 

P 

■ 

p 

c 

" 

■ ■ ! 

III 

ip 

Ipt 

lasaaaaoaaa 

• 

i 

IE 

: m i 

MM>1 

I 

■j 

P 

s 

1 

1 

| 

2 

■i 

§ 

; " i 

■ 

TTfr 

L*ll  1-1*1.’ 

i 

isaaaaaaaaa 

38309939983 

laaoaaaaaaa 

1 

| 

■ 

i 

i 

* 

! 

1 

:r:r 

Uli 

•W  r 

I i i [ vl  -i™nTimnin.33333333U.31.vl3333333331.LM.333U.L3333H.33030333333000 
i n.i  i lluminniiiminntTil  U-LLLLLU-s  ' ' LU.vU-333330303q0Q 33 30030 

iiiii  33M3BaaniD3tt333L3i  i vii  1 1 1 ii  inmi  i - i \i.nLLL3i]33]gjn333Ljmno 

NV  sU3333a33LL333LLLL'L3LLL33L3333LI.LL'.L'LLl.LLL333mJOaaBOOOaODOOaa3 
vs  LLL3LLL333U.3331  LL'  LLLsLLl.sLLU.LLss  sns  w\LL33D00030333a3333COQ 
ss.sssLLLLLsL3L3LL3L3LsLLL33LLL33LLLLLLLss  L'LLLILL33O009B0BDOOOOOOOEQD 
ssss\LL33lLLJ33L33333LLLL3LLLL3LsLsssLLsssLssLsLLL33DBaUUUU3aaSUOa3QBB 
sLssssL33LL3U3n333333LL3333LLJU333LLLLLLLsLL»LssLs33HBn3tl0a  30000000030 
r i xwsi  Tmnnnminma3333333L33a33333333vLsLLLLsLs^L3aa000003r,qaBEaaEBQ 
ll.\\>U3  33311BBBHHlinHIU33333LLL33LLLLLLLsssLLsLxLss330aBaBnUma:*;Ha0QBQ 
LL  • m i inmmiTiiiininmni  Lnnmi333LL\L\LL\LLLL\]jn«M  n iinnmm 


M 


LL3L'3333LLLL  3LLLLL3'L  \ s OOL3L"  • 3L' 3L3' vv  '3'L33  30^ 3L33vLLL0'LL 
'3LnLLLLLL3L'LLLLLLLLLLLL3LLLLLLLLL'3L"L'LL'L'LLL'LL3LLLLL'L'LLLL33L3 
LLL3LLL3LL  LLL'LLL3LL3LLLLLLLLLLLLLLLLL'  LL3LLLLLLLLLLL3V'LLLLL'L 3LLLLL 
LLLLLLLl.LL3LLLLL'L3LLLLLLLLLLL'LL3LLLL'  LLLLLLLLLL' LL333LL3L3L3LLLL3L3 
1L31LLLLLL331L3LLLLL3LLLL3LLLLLL  " LL3LL'LLLLLLLLLLLL33LL3L"'  L3LL  ~ 3L1L 
UI.LLLLJJLt  LLL 3LLL33LLLLLLL3LL'  L'LLL'LL'  L'LLLL'L3LLLLL'LLLLLL33LLLLLLL 
'3L'LL3LLLLL'L3L3LLLLLLLLLLLL  L LLL'  LLLLL"L'  L'LLLLL'L'L'LLLLLL3LLLL  ■ L 
JLLLL3LLLLLLLLL3LLLLLJLLLL'3'LLLLL"'LLLL'LLLLLLLLL'LL'LL'L33LLLLL  LLL 
3L"L3LLL3LL3L33LLL3LLILLLL3LLLLL'LLL3L'L'LLLL'LL'L33LLLL'LL3LLLLL3LL3 
LL3LL3LLLLLLLL3LLLL'LL'LLLLLLLLLLLLLLLLLLLL'L'LLL'L3'LLLLLLLLLL'LLL1LL 
3'LLLLLLLLLLLLL3L3'L'LLLLLL3L3LL'LLL'LLL'L'LLLLLLL3LL3LLLLLLL3LL33LLL3 
3LLLL-  'LLL3LL3JLLLL'LLLLLLLL3LLLL'L'L"L'LL"LLLLLnLLL'3L'LL'3LLLLLLLL 
'LLLLLLL JLLLLLL33LLLLLL'LLL'LLLLLLLLLL'3'LLLLLL3LLLLLLLLL3L'L3L'  L-LL3L 
3LL3LLLLLLLLLLL3LLLLL3LL3L'LLLL3L3L'LLLLLLLLLLL3LLLLLLL'LLLLLLLLLLL  LL 
3LLJLLLLJLLI.:  LL3L LLLLLL JLLLLLLLLLLL3LL3LLLLLLLLL'3LLLLLL'L3LLLLLL'3LLL 
LL1LL3LLLLL33LLL'LL'LLLL33L3'LLLLLL'LLL33L'LLLL3LLLLLLLL3LLL'3LLLLL'I.1 
1LLLL  '-LLLLLLLLLLL3LLL'LL3LLL3LLLLLLLLLLLLLL3'LLLLLLLLLLLLLLLLLLLLL3L 
n I.L  V L U '-LLLL3L3LLLL3LLL'3'3LL3LL3L3L3L33LLLLL3LL'3LLLL3LLL'LLL33L'L' 
V.LLL3L  -LLL'L3LLL3L3LLLLLLLLL3LLL3LL33'3LLL'LL'LLLLLLL33LL3LLLLLl  '..U. 
- 3'u  J'  I u :L1LJ1LLLLL3L3LL3L3333L3L3''LLL3LLLLLLLLLLLL'3LLL33L3LL • LL  3L 

I '■  1 y-  • JLL1L3'LLLL3LLLLL3LLLLL3U1L333LL3LLLLL'L3LLL3LLLL33'L'LLLL 

J...LLL  JLL33L33LLLLL3333L33L3LL3LLLLLLLLLLLLL3'L3LLLL3'LLLLLLLLLL 
1 ’ . ..  II  : L)LLLLL3J33LL3La3333L330333LLL3L'33L33LLL3'LLLL3LLJLL'LL'.L 

II  V..  I ML  11  3L33333LLL3L1M33LLLL33'L3L3LL3L3L3LLL3LLLL3L'33LLLLL  LL  » 
'-Uu  II.  Jl.  LL  J3LLLLL33LLLL3MM33LL33L3LLLLLL333L3LLLLLL3LLLLLLL3LLLL'  t 1 
1L-  J . 3t  L:.L3LLLLLLL'3'L3BMMBLLLLLLLL3L3LLJ3LLLL3L3LLLLLLLLLJ3LLLLLL 

J It  1.1'  LL  JLl.JLL3LLL33LLL3SMM333LL3LLLLLLL3LLLL3L3LLL3LLLLL3LL'3LLL'  LL 
-y.  1 LL!.:  11.  1LLL3L3LL3LL3LLLM«3'L3'33L3LLLLL3LLJ3LL3L'LLL3LLLLL3LLL  LLLL 
L L _LLLJL.11L313LL3L3'3LLLL31D'33LL3L3LLL3L3LLJ33LLL3'3LLLL33LLLLLL3  1L1 
LLL3L3LL3LLLLLL33LLL3LL390QL3L3L'3L3LLL3LLLL33L33LLL3L33LL3'3LL  'LLL1.L 
' L3LL3L3L3LLLLL333L3LL33La0'3LL3L3L3333L'3L3LLLL3LLL'L3LLL3LLL33LLL31 J 
njLLL'.LLLL3L3333  0L3L3LL3LL3LL3LLLLL3LLLLL3LLLLL3L3L3LLL3333LLL3'L1..J 
-JJ  - JI.L  3 31L33L3L33LL3QL3LLL33LL3L33L3L3'LL3L3LLL3LLL3LLL3L3L3LJL  JL 
3L31L3LL33333L3LL333LLL3HJLLLLLLLL33LL3L'L3LL33LL33LLL'LLL333LLLLLLLLL 


Figure  8.16  Frame  6 after  the  hybrid  filter 
(TBRQ  = 15  dB) 
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Figure  8.17  Frame  10  after  the  hybrid  filter 
(TBRq  = 15  dB) 
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(example  3) 
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This  example  is  similar  to  the  previous  one,  but  the 
hybrid  instead  of  the  recursive  filter  is  used.  The  frame 
1 is  not  shown,  since  both  filters  reduce  to  a 2D-recursive 
filter  for  this  frame.  Comparing  the  results  we  can  see 
that  the  recursive  filter  (example  1)  presented  better 
performance  than  the  hybrid  one.  As  expected,  the  hybrid 
filter  had  almost  constant  performance  (15  dB  of  processing 
gain)  for  all  frames,  since  it  has  no  transient  in  time. 
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Example  3 

X / 

Filter:  Recursive 


TBR^  = -3  dB 

» Threshold  = 0.35 

Figures:  8.18  through  8.23 

1 

-•  ' The  filter  performance  for  other  frames  is  shown 

in  table  8.4,  where  the  effect  of  the  transition  in  time 
at  the  first  5 frames  is  clearly  shown. 

Example  4 

Filter:  Hybrid 
TBRi  = -3  dB 

Threshold  = 0.4 

Figures:  8.24  and  8.25 


This  example  is  similar  to  the  previous  one,  but 
using  the  hybrid  instead  of  the  recursive  filter.  The 
processing  gain  was  around  12.5  dB  for  all  frames.  From 
figures  8.24  and  8.25  we  can  see  that  TBR . = -3  dB  is  about 
the  limit  of  detection  for  the  hybrid  filter. 

4 . Target  Tracking 

In  what  follows  we  present  some  examples  of  target 
J detection  and  tracking  from  frame  to  frame.  The  background 

has  the  same  characteristics  as  before.  The  first  five 
frames  are  entirely  processed  in  order  to  acquire  the  target 
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Figure  8.18  Original  image  at  frame  1 (TBI^  = -3  dB) 
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Figure  8.19  Frame  1 after  the  recursive  filter 
(TBRq  = 9 dB) 
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Figure  8.20  Original  image  at  frame  6 (TBR.  = -3  dB) 
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Figure  8.21  Frame  6 after  the  recursive  filter 
(TBRQ  =15.8  dB) 
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Figure  8.22  Original  image  at  frame  10 
(TBRi  = -3  dB) 
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Figure  8.23  Frame  10  after  the  recursive  filter 
(TBR0  =15.9  dB) 
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Figure  8.24  Frame  6 after  the  hybrid  filter 
(TBRq  =12.5  dB) 
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Figure  8.25  Frame  10  after  the  hybrid  filter 
(TBRq  =12.5  dB) 
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(search  mode) . The  window  mode  is  implemented  at  frame  6 
and,  thereafter,  the  pictures  are  processed  only  in  a 
small  20  x 20  window  centered  at  the  target  predicted 

position.  In  the  examples  we  have  created  binary  pictures 

' 

to  emphasize  the  performance  of  the  filters. 

Example  5 

Filter:  Recursive 

TBRi  = 0 dB 

Threshold  = 0.5 

Figures:  8.26  through  8.32 

Example  6 

Filter:  Hybrid 

TBRi  = 0 dB 
Threshold  = 0.5 
Figures  8.33  through  8.36 

Observe  that  the  performance  of  the  hybrid  filter 
is  inferior  to  that  of  the  recursive  filter  in  example  5. 

It  can  also  be  seen,  in  figure  8.36,  that  the  tracking  is 
better  after  the  first  5 frames.  The  reason  is  that  in  the 
window  mode  the  measurement  error  of  the  centroid  position 
is  limited  by  the  size  of  the  window. 
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Figure  8.26  Original  image  at  frame  6 


(TBR^  = 0 dB) 
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Figure  8.27  Frame  6 after  the  recursive  filter 
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Figure  8.28  Original  image  at  frame  12 
( TBRi  = 0 dB) 
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Figure  8.29  Frame  12  after  the  recursive  filter 
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Figure  8.30  Original  image  at  frame  18 
(TBRi  » 0 dB) 
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Figure  8.31  Frame  18  after  the  recursive  filter 
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Figure  8.32  'Target  tracking  with  the 

recursive  filter  (example  5) 
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Figure  8.33 


Frame  6 after  the  hybrid  filter 
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Figure  8.34  Frame  12  after  the  hybrid  filter 
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Figure  8.35  Frame  18  after  the  hybrid  filter 
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Example  7 


Filter:  Recursive 


TBRi  = -3  dB 


Threshold  = 0.4 


Figures:  8.37  through  8.43 


Example  8 


Filter:  Hybrid 


TBRi  = -3  dB 


Threshold  = 0.4 


Figures:  8.44  through  8.47 

Comparing  examples  7 and  8 we  see  that  both  filters 
were  able  to  track  the  target,  but  the  recursive  filter  is 
clearly  superior  to  the  hybrid  one. 
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Figvire  8.37  Original  image  at  frame  6 (TBR^  = -3  dB) 
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Figure  8.38  Frame  6 after  the  recursive  filter 
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Figure  8.39  Original  image  at  frame  12  (TBR^  = -3 
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Figure  8.40  Frame  12  after  the  recursive  filter 
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Figure  8.41  Original  image  at  frame 
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Figure  8.42  Frame  18  after  the  recursive  filter 
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Figure  8.43  Target  tracking  with  the 

recursive  filter  (example  7) 
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Figure  8.44  Frame  6 after  the  hybrid  filter 
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Figure  8.45  Frame  12  after  the  hybrid  filter 
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Figure  8.46  Frame  18  after  the  hybrid  filter 
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IX.  CONCLUSIONS 
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A.  SUMMARY  OF  RESULTS 

The  experiments  with  several  real  life  pictures  allow 
us  to  draw  two  important  conclusions  concerning  image 
modeling.  First,,  the  hypothesis  of  an  autocorrelation 
function  with  separable  kernels  in  the  vertical  and  hori- 
zontal directions  worked  quite  well.  This  conclusion  is 
very  important,  because  the  derivation  of  the  dynamic  model 
is  substantially  simplified  by  avoiding  the  difficult 
problem  of  two-dimensional  spectral  factorization.  Second, 
the  first-order  model,  used  by  most  researchers,  is  a poor 
approximation  for  pictures  with  few  details,  although  it  is 
a very  good  model  for  pictures  with  many  details.  The  auto- 
correlation function,  suggested  in  this  research,  seems  to 
be  a good  choice,  because  it  includes  the  first  and  second 
order  models,  as  particular  cases,  and  also  permits  a simple 
method  of  parameter  identification  which  can  be  easily 
implemented. 

The  performance  evaluation  of  the  sub-optimum  two- 
dimensional  recursive  filters  of  Habibi’s  type  [3-5]  allow 
us  to  draw  some  interesting  conclusions.  Although  sub- 
optimum, these  filters  are  not  so  far  from  optimality.  In 
the  worst  case,  the  error  variance  was  around  15%  greater 
than  the  one  obtained  by  the  optimum  interpolator. 
Rosenfeld's  filter  [3]  is  the  best  and  presented  degradation 
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around  8%.  Shachar's  filter  [5]  has  the  simplest  algorithm 
for  gain  calculation  which  uses  a very  good  approximation, 
since  its  performance  was  quite  close  to  that  of  a similar 
filter  suggested  here  which  calculates  the  gains  without 
approximation - 

The  reason  why  Rosenfeld's  filter  [3]  is  the  best  is 
because  the  estimation  error  is  orthogonal  to  the  obser- 
vation of  the  pixel  under  estimation,  although  this  does 
not  happen  for  the  rest  of  the  data  set.  This  conclusion 
leads  directly  to  the  hybrid  filters  introduced  in  this 
research.  From  the  results  shown  in  Table  5.1,  it  can  be 
seen  that  the  non-recursive  filter  [6] , suggested  by  Bar- 
Yehoshua  [6],  although  using  only  9 observations,  presents 
better  performance  than  the  optimum  (in  the  sense  of 
recursive  filtering)  interpolator,  for  correlations  as  high 
as  0.94.  The  conclusion  is  that  most  of  the  information 


about  the  pixel  gray  level  under  estimation  resides  on  its 
closest  neighbors.  For  this  reason,  the  hybrid  filter 
presents  the  best  performance,  because  it  makes  use  of  all 
the  observations  in  the  neighborhood  and  also  those  farther 
observations  used  by  the  recursive  filter.  Experiments 
with  real  life  pictures  (figures  5.12  and  5.13)  show  the 
ability  of  the  hybrid  filter  in  smoothing  out  observation 
noise  in  pictures  with  signal-to-noise  ratio  as  low  as 
-3  dB.  The  processing  gain  for  this  case  was  11.1  dB,  which 
compares  favorably  with  previous  results. 
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For  images  sequenced  in  time,  the  three-dimensional 
recursive  and  hybrid  filters  were  developed  in  order  to 
exploit  the  correlation  in  time.  The  recursive  filter  is 
adequate  for  the  cases  where  the  scene  is  the  same  during 
several  frames,  otherwise  the  exploitation  of  the  time 
correlation  is  not  advantageous,  due  to  the  filter  transient 
at  the  first  few  frames.  Its  structure  is  intuitively 
appealing,  since  it  reduces  to  the  one-dimensional  Kalman 
filter  for  the  case  of  no  spatial  correlation,  and  to  a two- 
dimensional  recursive  filter  for  the  case  of  no  correlation 
in  time.  Numerical  results  show  that  the  use  of  time  corre- 
lation can  be  quite  advantageous  for  the  case  of  time 
correlation  greater  than  spatial  correlation  (error  variance 
is  reduced  by  33%) , but  it  is  not  so  advantageous  for  the 
opposite  case  (error  variance  is  reduced  by  6%) . The  hybrid 
filter  is  adequate  for  the  case  where  the  scene  changes  at 
every  two  frames,  since  it  makes  use  of  only  the  previous 
and  the  present  frames.  Of  course,  it  can  also  be  designed 
for  the  case  of  a fixed  scenario  (see  Chapter  V)  in  which 
case  it  improves  the  performance  of  the  recursive  filter. 

To  detect  targets  from  cluttered  background  images  a 
likelihood  ratio  was  constructed  which  uses  the  observation 
of  the  pixel  and  the  background  prediction  for  the  pixel 
gray  level.  A quite  general  situation  was  considered. 

The  target  and  background  are  considered  as  two  kinds  of 
textures.  The  image  is  modeled  as  a weighted  combination  of 
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three  additive  components:  target,  background  and  obser- 
vation  noise.  Therefore,  it  can  accommodate  the  case  of 
infrared  pictures  and  also  visible  pictures,  where  there  is 
a replacement  process,  target  or  background.  The  prediction 
of  the  background  is  given  by  the  recursive  and  hybrid  filter 
or,  also,  by  a non-recursive  filter.  The  likelihood  ratio 
was  particularized  for  four  special  cases  of  interest,  which 
had  their  performance  analyzed.  The  likelihood  ratio  enhanced 
the  target- to-background  ratio.  The  target  is  extracted  by 
threshold  detection.  The  threshold  is  chosen  in  order  to 
maintain  the  probability  of  false  alarm  below  some  desired 
value.  The  performance  of  the  detection  process  is  highly 
related  to  the  background-to-observation  noise  ratio,  since 
this  ratio  is  responsible  for  the  background  prediction. 

The  target  detection  and  tracking  problem  was  simulated 
using  comptuer  generated  images  with  characteristics  similar 
to  those  of  real  life  images  used  in  the  other  experiments. 
Both  three-dimensional  recursive  and  hybrid  predictors  were 
used  and  compared.  The  likelihood  ratio,  using  both  filters, 
was  able  to  detect  a target  at  target- to-background  ratios 
as  low  as  -3  dB.  The  processing  gain  was  around  19  dB,  for 
the  recursive  filter,  and  around  12.5  dB  for  the  hybrid 
filter.  The  recursive  filter  presented  better  performance, 
as  should  be  expected  from  the  results  of  Tables  5.3  through 
5.5,  where  its  prediction  is  better  than  that  of  the  hybrid 
filter,  for  the  case  of  high  background-to-noise  ratio. 
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B.  SUGGESTIONS  FOR  FUTURE  RESEARCH 

Many  pictures  are  clearly  non-homogeneous , therefore 
the  development  of  space-varying  models  is  the  natural 
way  to  proceed  in  order  to  improve  the  model  suggested  in 
this  research.  The  recursive  techniques  are  directly 
amenable  to  the  analysis  of  space-varying  models.  The 
reliance  on  "a  priori"  information  is  a very  important 
area.  The  method  of  parameter  identification  suggested 
here  is  applicable  for  the  case  of  constant  parameters  for 
the  whole  frame.  In  the  non-homogeneous  case  such  identi- 
fication has  to  be  accomplished  pixel  by  pixel,  in  order 
to  improve  the  robustness  of  this  statistical  approach 
in  face  of  modeling  errors. 

Due  to  the  enormous  computational  load  of  the  optimum 
multi-dimensional  recursive  filters,  and  also  because  the 
observations  far  away  carry  negligible  information  about 
the  pixel  to  be  estimated,  the  natural  way  is  to  look  for 
sub-optimal  recursive  filters  that  require  less  computation. 

The  hybrid  filters  suggested  here  can  be  equally  applied 
with  other  recursive  filters,  as,  for  example,  the  reduced 
dimension  filters  proposed  by  Panda  and  Kak  [9]  and  Woods 
and  Radiwan  [11] . 

The  decision  rule  suggested  here  is  directed  to  the 
detection  of  targets  pixel  by  pixel,  in  order  to  be  indepen- 
dent of  the  target  shape.  The  choice  of  the  threshold  was 
based  on  the  probability  of  a false  alarm  of  each  isolated 
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